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Executive  Summary 


Reliable  transmission  of  information  is  an  important  issue  in  military  and  com¬ 
mercial  communication  systems.  Even  though  modern  techniques  provide  some 
error  performance  improvement,  systems  experiencing  large  interference  can  un¬ 
dergo  unacceptable  degradation.  The  principal  investigator  and  his  research 
team,  under  previous  contracts  with  the  United  States  Air  Force  through  Rome 
Laboratory,  developed  a  system  that  significantly  improves  performance  in  sev¬ 
erely  degraded  channels.  This  report  represents  the  culmination  of  a  study 
focusing  on  more  extensive  characterizations  of  a  robust  autoregressive  locally 
optimum  (ARLO)  detector  operating  in  a  variety  of  noise  environments. 

In  a  communication  system,  undesired  signals  (such  as  channel  noise,  jam¬ 
mers,  or  cross-channel  interference)  are  a  primary  cause  of  signal  corruption 
and  performance  degradation.  If  the  channel  were  corrupted  with  only  additive 
white  Gaussian  noise,  the  chosen  detector  would  be  a  linear  correlator,  as  it  is 
an  optimum  detector  for  this  scenario.  However,  a  linear  correlator  is  usually 
sub-optimum  for  other  types  of  noise,  such  as  would  be  encountered  in  an  inten¬ 
tionally  or  unintentionally  hostile  environment  or  in  a  multi-user  application.  In 
this  scenario,  traditional  techniques  usually  precede  the  linear  correlator  with 
a  fixed  nonlinearity  [1],  whose  purpose  is  to  negate  the  effects  of  the  additional 
noise,  and  thus  improve  the  performance  of  the  linear  correlator.  The  ARLO 
algorithm,  in  contrast,  uses  an  adaptive  nonlinearity  that  utilizes  probability 
density  function  (pdf)  estimation  techniques  to  alleviate  the  need  for  a  priori 
knowledge  of  the  channel  noise  characteristics.  The  mathematical  basis  for  the 
algorithm  is  as  follows. 

Consider  a  system  that  has  a  received  signal  vector,  r,  given  by  the  following 
expression, 

r  =  STO-|-n,  (1) 


where  s™  is  the  transmitted  signal,  m  =  Oorl,  and  n  is  the  total  noise  signal 
(interference  plus  background  noise). 

Letting  the  observed  value  of  r  be  p,  the  LO  nonlinearity  used  in  the  ARLO 
detector  is. 


9i  (p) 


/n(/3) 


(2) 


where  /„  (•)  is  the  joint  pdf  of  the  noise. 
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Therefore,  the  LO  detector  likelihood  function,  I  {p) ,  for  the  detection  of 
known  binary  signals  can  be  written  as, 

choose  Hi 

lip)  =  j:li[isii-soi)gi{p)]  <  7  (3) 

choose  Ho 

where  7  is  an  appropriately  chosen  decision  threshold. 
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Chapter  1 

Introduction 


In  modern  military  and  commercial  communication  systems,  it  is  imperative 
that  information  be  transferred  in  a  reliable  and  secure  manner.  Current  tech¬ 
niques,  such  as  those  used  in  spread  spectrum  (SS)  systems,  provide  a  measure 
of  security  and  reliability.  Systems  that  undergo  extreme  interference,  however, 
can  exhibit  unacceptable  performance.  It  is,  therefore,  imperative  to  investi¬ 
gate  techniques  to  combat  the  degradation  experienced  by  such  systems.  Since 
the  decrease  in  acceptability  is  caused  by  interference,  it  is  logical  to  develop 
methods  to  mitigate  these  undesirable  signals. 

Locally  optimum  (LO)  detection  is  a  technique  that  lends  itself  to  such  sit¬ 
uations.  Its  performance  is  close  to  that  of  the  optimum  Bayesian  or  globally 
optimum  (GO)  detector  with  a  much  reduced  complexity.  It  is  also  robust  in 
that  it  can  be  implemented  to  adapt  to  changing  environments.  Previous  work 
[1]  has  focused  on  such  LO  detectors  that  utilize  estimates  of  the  probability 
density  function  (pdf)  alleviating  the  need  for  any  a  priori  knowledge  of  the 
noise  characteristics.  Multidimensional  signal  processing  can  improve  perfor¬ 
mance  in  these  systems.  However,  the  complexity  increases  enormously  as  the 
dimensionality  of  the  joint  pdf  (jpdf)  of  the  signaling  environment  grows.  More 
recently,  techniques  have  centered  on  the  use  of  autoregressive  (AR)  modeling 
methods  to  better  estimate  the  noise  environment  and  enhance  performance. 
The  Autoregressive  Locally  Optimum  (ARLO)  detector  is  one  that  successfully 
combines  these  techniques.  The  use  of  AR  process  models  greatly  reduces  the 
dimensionality,  and  hence  the  complexity,  of  the  resulting  ARLO  detector.  The 
statistics  of  the  jpdf  remain  an  issue.  To  this  end,  pdf  estimation  techniques 
for  independent  identically  distributed  (iid)  noise  samples  prove  useful  since  the 
input  sequence  to  the  AR  model  is  iid  albeit  characterized  by  some  unknown 
pdf.  The  necessary  AR  model  parameters  are  determined  by  well-known  spec¬ 
tral  estimation  techniques.  Thus,  the  goal  of  this  research  effort  is  to  further 
characterize  performance  improvements  in  robust  LO  detectors  using  AR  tech¬ 
niques. 

In  pursuit  of  this  goal,  simulations  were  initially  performed  using  MATLAB® 
and  later  re-coded  into  C.  The  simulation  uses  B  total  received  symbols  sampled 
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at  a  rate  of  Nb  samples/symbol.  Standard  AR  model  estimation  techniques  are 
applied  to  all  of  the  received  data.  The  bit  error  performance  (BER)  of  five 
detectors  are  computed  and  compared.  The  detectors  studied  are  the  linear 
correlator,  the  LO  detector  without  memory,  the  linear  correlator  and  the  LO 
detector  operating  on  whitened  noise,  and  the  ARLO  detector.  A  block  diagram 
illustrating  this  system  is  shown  in  Figure  1.1. 


Nonlinearity 


(a) 


pdf  Estimate  I 


(b) 


Figure  1.1:  (a)  ARLO  system  block  diagram,  (b)  Expanded  view  of  Estimate 
h’  block. 

The  interference  environment  was  generated  using  an  AR  process  by  passing 
iid  noise  through  a  given  AR  model  (i.e.,  the  {at}  are  known).  The  communi¬ 
cation  system  employed  binary  phase  shift  keyed  (BPSK)  signals  so  that  each 
symbol  represented  a  bit.  The  AR  process  model  used  a  Id***  order  filter.  The 
spread  spectrum  system  used  a  processing  gain  of  16  chips/bit  and  the  number 
of  samples/chip,  Nc,  was  set  to  4.  The  noise  pdf  was  estimated  using  the  his¬ 
togram  method.  A  Monte  Carlo  simulation  was  performed  using  a  number  of 
iterations. 

The  organization  of  this  report  is  as  follows.  First,  pdf  estimation  tech¬ 
niques  are  evaluated  for  performance  enhancement  in  Chapter  2.  Next,  the 
performance  of  the  ARLO  detector  is  compared  with  other  standard  detectors 
in  multiple  interference  environments  in  Chapter  3.  Additionally,  this  chapter 
reviews  the  effect  of  detector  variables  on  performance.  Chapter  4  focuses  on 
the  conversion  of  the  simulation  into  a  high  level  language.  In  addition,  this 
chapter  analyzes  the  execution  speed  of  the  simulation  and  reviews  the  effect  of 
detector  parameters  on  said  speed.  Finally,  a  summary  of  the  research  study 
is  provided  in  Chapter  5,  including  a  discussion  concerning  future  research  and 
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development  of  the  ARLO  detector. 
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Chapter  2 

Probability  Density  Function 
Estimation 


The  current  form  of  the  simulation  uses  the  histogram  method  for  estimating 
the  pdf  of  the  whitened  received  signal.  The  derivative  of  the  resulting  pdf  is 
then  used  to  generate  the  decision  function,  g  {p),  for  choosing  between  the  two 
BPSK  symbols.  Simulation  results  have  revealed  that  a  large  portion  of  the 
processing  time  is  devoted  to  the  calculation  of  g  {p).  Therefore,  the  goal  was  to 
determine  if  a  faster,  more  efficient  alternative  to  the  histogram  method  could 
be  found  for  estimating  the  pdf  of  the  received  signal. 

Efforts  became  focused  on  estimating  the  unknown  pdf  of  a  correlated  data 
sequence  via  POCS  methods.  However,  after  initial  research,  it  was  determined 
that  the  POCS  method  may  not  offer  the  coding  simplicity  and  speed  required 
due  to  its  potential  complexity.  Research  was  then  concentrated  on  using  arti¬ 
ficial  neural  networks  (ANN)  for  pdf  estimation  [5].  This  work  proposes  that  a 
mixture  of  radial  basis  functions  (RBF)  can  approximate  an  unknown  pdf.  An 
RBF  ANN  can  be  trained  in  such  a  way  that  it  outputs  the  approximated  pdf, 
even  when  the  input  data  is  correlated.  However,  in  this  form,  the  ANN  method 
was  more  complicated  than  necessary  due  to  the  fact  the  actual  simulation  data 
is  uncorrelated  because  the  received  signal  is  whitened  prior  to  the  determina¬ 
tion  of  its  pdf.  This  fact  simplifies  the  problem  at  hand,  thus  making  it  possible 
to  use  a  similar  but  less  complicated  theory  called  Parzen’s  Estimator  [4]. 

2.1  Histogram  Method 

The  histogram  method  attempts  to  estimate  the  pdf  of  a  signal  through  the 
use  of  a  histogram  to  approximate  its  continuous  pdf.  Each  signal  sample  is 
deposited  into  the  appropriate  bin,  based  on  its  value.  This  process  creates  a 
histogram  which  describes  the  pdf  of  the  signal.  The  greater  the  number  of  bins, 
the  better  the  histogram  represents  the  continuous  pdf  of  the  actual  received 
signal.  This  method  is  very  simple  to  implement,  and  quite  accurate  relative  to 
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its  level  of  complexity. 


2.2  Parzen’s  Estimator 

Parzen’s  estimator  attempts  to  estimate  a  pdf  through  the  use  of  radial  basis 
functions  (RBF).  RBFs  are  parametric  functions  which  satisfy  the  axiomatic 
properties  of  pdfs.  They  assume  a  maximum  value  at  zero  radius  and  mono- 
tonically  decay  from  this  maximum  as  the  radius  increases  [5].  For  instance, 
the  Gaussian  density  function  is  considered  to  be  an  RBF.  Parzen’s  estimator 
proposes  that  the  estimated  pdf  of  a  set  of  correlated  data  is  the  weighted  sum 
of  individual  RBFs  centered  at  each  data  point.  This  can  be  written  as, 

/x  ,  (2.1) 

where  n  is  the  number  of  data  points,  k  (•)  is  the  RBF,  and  h  (n)  is  a  smoothing 
factor  [4].  Currently,  the  Gaussian  function  is  being  used  as  the  RBF,  fc(-), 

and  h  (n)  is  the  inverse  root  of  the  sample  size  [4], 

h  (n)  =  n~ .  (2.3) 

The  above  definition  of  the  smoothing  factor  results  in  the  fact  that  as  the 
value  of  p  increases,  the  bias  of  the  estimator  increases,  while  the  variance  of 
the  estimator  decreases.  If  p  takes  on  a  value  of  2,  the  resulting  estimator  is 
unbiased.  As  stated  above,  the  derivative  of  the  estimated  pdf  must  be  taken  in 
order  to  determine  the  g  {p)  function,  so  achieving  a  smooth  pdf  is  important. 
However,  the  resulting  unbiased  estimator,  when  p  =  2,  is  not  smooth.  Smooth¬ 
ness  is  achieved  when  the  value  of  p  is  greater  than  3,  but  these  estimators 
are  biased.  Comparatively,  in  the  histogram  method,  the  number  of  bins  used 
affects  the  smoothness  of  the  resulting  pdf.  Currently,  the  simulation  employs 
a  histogram  with  thirty-three  bins.  This  results  in  a  jagged  pdf.  Smoother 
results  are  achieved  when  the  number  of  bins  approaches  twelve.  The  various 
values  of  p  yield  approximately  the  same  estimate  of  the  pdf.  When  p  =  3.4, 
the  lowest  MSE  is  achieved.  However  it  is  only  slightly  smaller  than  the  other 
MSE  values.  All  of  the  histogram  estimates  are  somewhat  jagged  at  transition 
points,  thereby  resulting  in  large  MSEs. 

The  Parzen  estimator  can  prove  to  be  a  reliable  estimator  of  the  pdf  of  the 
whitened  received  signal  data.  The  method  results  in  smooth,  estimated  pdfs 
with  low  MSEs.  In  order  to  make  a  definitive  decision  on  whether  the  Parzen 
estimator  is  more  effective  than  the  histogram  method,  further  research  needs 
to  be  conducted.  The  issues  which  need  to  be  addressed  are  the  following: 
the  exact  estimation  time  of  each  method,  the  capabilities  of  each  method  to 
reliably  estimate  pdfs  using  fewer  data  points,  and  the  capacity  of  each  method 
to  estimate  other  types  of  non-Gaussian  pdfs. 
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Chapter  3 

Performance  Analysis 


The  ultimate  measure  of  the  total  system  performance  of  a  communication  sys¬ 
tem  is  the  BER,  the  ratio  of  erroneous  bits  received  to  the  total  number  of  bits 
transmitted.  In  order  to  be  able  to  measure  the  BER,  a  SS  communication 
system  was  previously  implemented  using  MATLAB®  as  a  simulation  tool  [3]. 
The  simulation  will  allow  us  to  study  the  BER  performance  of  the  system  in 
varying  interference  environments  as  well  as  under  different  configurations  of 
system  parameters.  The  first  section  of  this  chapter  will  describe  the  various 
interference  environments  implemented  in  the  simulation.  The  second  section 
will  analyze  the  detector  performance  against  each  interferer.  The  ARLO  detec¬ 
tor  has  a  number  of  parameters  whose  values  can  affect  the  performance  of  the 
communication  system.  Therefore,  the  third  section  of  this  chapter  will  discuss 
in  detail  the  effects  of  choosing  alternate  parameter  values. 

3.1  Interference  Environments 

Eour  different  interferers  were  implemented  for  use  in  the  simulation,  continu¬ 
ous  wave  (CW),  partial  band  (PB),  dual  continuous  wave  (2CW),  and  mixed 
interferers. 

3.1.1  Continuous  Wave 

The  CW  interferer  is  simply  a  sine  wave  of  a  predetermined  frequency,  with 
zero  phase  shift,  as  shown  in  Eigure  3.1.  The  amplitude  of  the  sine  wave  is 
determined  by  the  interference  to  signal  ratio  (ISR)  as. 


Amplitude  = 


(3.1) 


where  I  is  the  number  of  CW  interferers  (one,  in  the  single  CW  case). 
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Masnitude 


Figure  3.1:  CW  Interferer,  100  data  samples. 


The  frequency  of  the  sine  wave  is  determined  as, 

/  =  NeB,  (3.2) 

where  Nc  is  the  number  of  samples  per  chip  and  B  is  the  number  of  symbols 
transmitted  over  the  channel.  The  power  spectra  of  a  CW  interferer  with  ISR  = 
30dBW  is  shown  in  Figure  3.2. 

3.1.2  Partial  Band 

The  PB  interferer  exhibits  a  large  degree  of  power  over  a  chosen  band  of  frequen¬ 
cies.  This  jammer  is  created  by  passing  an  independent,  identically  distributed 
(IID)  white  noise  process  through  a  predetermined  AR  filter,  which  is  created 
by  taking  the  polynomial, 

poly  =  [1  -  0.75] ,  (3.3) 
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Power  Spectra  for  CW  Interference  Environment 


Frequency  {with  respect  to  sampling  frequency) 


Figure  3.2:  CW  power  spectra. 


and  convolving  it  with  itself  four  times.  This  “auto-convolution”  process  creates 
a  16*^-order  filter  with  filter  coefficients, 


-1 

12 

-67.5 

236.25 

-575.859 

1036.55 

-1425.25 

1527.06 

-1288.45 

858.969 

-450.959 

184.483 

-57.651 

13.3041 

-2.13815 

0.213815 

-0.0100226 


(3.4) 


The  power  spectra  of  the  filter  described  by  (3.4),  with  ISR  =  30  dBW,  is  shown 
in  Figure  3.3. 
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Power  Spectra  for  PB  Interference  Environment 


Figure  3.3:  PB  power  spectra. 


3.1.3  Dual  Continuous  Wave 

The  2CW  interferer  is  easily  created  by  summing  two  single  CW  interferers. 
The  amplitude  of  each  sine  wave  is  determined  as  shown  in  (3.1),  with  1  =  2. 
The  frequency  of  each  sine  wave  is  calculated  by, 


NeB 

0.3NcB  ’ 


(3.5) 


with  the  variables  the  same  as  in  (3.2).  Figure  3.4  shows  the  2CW  interferer, 
with  ISR  =  30dBW,  along  with  its  two  component  sine  waves.  The  power 
spectra  for  the  2CW  interferer  is  shown  in  Figure  3.5. 


3.1.4  Mixed 

The  mixed  interferer  is  created  by  summing  a  CW  interferer  with  a  PB  interferer. 
The  ISR  of  the  combination  is  fixed;  the  sum  of  the  ISR’s  of  the  components 
equals  the  value  chosen.  The  power  spectra  of  such  a  jammer  is  shown  in  Figure 

3.6. 
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Figure  3.4:  2CW  interferer,  shown  with  its  component  parts. 


Power  Spectra  for  2CW  Interference  Environment 


Figure  3.5:  2CW  power  spectra. 
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Power  Spectra  for  MX  Interference  Environment 


Figure  3.6:  MX  power  spectra. 


3.2  Detector  Performance 

The  SS  communication  system  simulation  has  a  number  of  parameters  that  de¬ 
fine  the  system  operation.  In  order  to  properly  analyze  the  performance  of  the 
ARLO  detector,  and  compare  it  to  the  BER  performance  of  other  detectors, 
default  values  for  each  parameter  must  be  chosen.  Assessing  the  effect  on  per¬ 
formance  of  changing  these  parameters  is  the  focus  of  the  third  section  of  this 
chapter.  Default  values  were  chosen  based  on  various  factors  such  as  execution 
time  and  typical  values  encountered  in  practice.  Table  3.1  shows  the  default 
variables  chosen  for  simulation  purposes.  All  simulations  can  be  assumed  to 
have  been  executed  with  these  parameter  values,  unless  otherwise  noted. 

Four  other  detectors  are  implemented  in  the  simulation  for  comparison  with 
the  ARLO  detector.  The  simplest  is  the  linear  correlator  (LC),  which  is  the 
optimum  detector  in  an  additive  white  Gaussian  noise  environment,  such  as  the 
thermal  background  noise  used  in  the  simulation.  The  next  detector  is  the  LC 
on  Whitened  Noise  (LCW).  This  detector  performs  as  a  linear  correlator,  but  it 
operates  on  the  received  signal  data  stream  after  it  has  been  whitened  to  reduce 
the  noise  component.  The  third  detector  is  the  LO  detector,  which  is  the  basis 
for  the  ARLO  detector.  The  final  detector  is  the  LO  on  Whitened  Noise  (LOW). 
This  detector  performs  as  a  LO  detector,  but  it  operates  on  the  received  signal 
data  stream  after  it  has  been  whitened  to  reduce  the  noise  component.  The 
following  sections  will  analyze  the  performance  of  the  ARLO  detector  within 
each  interference  environment,  including  a  scenario  with  no  interferer  present, 
using  the  default  system  parameters  as  shown  in  Table  3.1. 
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Type 

Parameter 

Parameter  Description 

Value 

pdf 

Estimation 

Bh 

Total  symbols  to  calculate  histogram 

128 

K 

Data  points  to  calculate 

33 

T 

Range  of  support  for  pdf 

1000 

Spread- 

Spectrum 

Nc 

Samples  per  chip 

4 

PG 

Processing  Gain  (chips  per  bit) 

16 

Nb 

Samples  per  symbol 

(PC'  *  Nc) 

B 

Symbols  transmitted/received 

16,384 

Ntot 

Total  samples  in  data  stream 

mummisusm 

AR-model 

Bar 

Total  symbols  to  estimate  AR  parameters 

MTTm 

Bg 

Total  symbols  to  calculate  g  (.) 

Bh 

P 

Number  of  AR  coefficients 

16 

Parameters 

sigmat2 

Variance  of  thermal  noise 

1 

limit 

Limiter  on  noise 

hard 

Interferer 

ISR 

Interferer-to-signal  ratio 

30  dB 

Table  3.1:  Default  parameters  for  communication  system  simulations.  All  simu¬ 
lations  can  be  assumed  to  have  been  executed  with  these  parameter  values,  unless 
otherwise  noted. 

3.2.1  Partial  Band  Interferer 

The  PB  jammer  provides  for  an  extremely  hostile  environment  for  a  communi¬ 
cation  system  channel,  as  shown  in  the  power  spectra  plot  in  Figure  3.8.  The 
interferer  overpowers  the  signal  by  nearly  40  dBW  at  peak  jammer  power,  and 
further  overloads  the  signal  over  a  significant  portion  of  the  signal’s  main  lobe. 
As  can  be  seen  from  Figure  3.7,  the  ARLO  detector  greatly  outperforms  the 
other  detectors,  even  the  LO  detector  upon  which  it  is  based.  At  a  signal-to- 
thermal-noise  ratio  (SNR)  of  0  dB,  the  ARLO  detector  is  able  to  perform  at  a 
BER  of  less  than  10^^,  or  less  than  1  bit  error  for  every  100  transmitted,  despite 
the  powerful  PB  jammer.  The  next  best  detector  returns  a  BER  of  greater  than 
10^^,  or  approximately  19  bit  errors  for  every  100  transmitted.  In  this  scenario, 
ARLO  is  clearly  the  best  detector  of  all  those  tested. 

3.2.2  Continuous  Wave  Interferer 

In  the  CW  interference  environment,  it  can  be  seen  from  Eigure  3.9  that  the 
ARLO  detector  performs  adequately,  albeit  slightly  worse,  than  the  other  de¬ 
tectors  in  the  less-hostile  CW  environment,  with  the  exception  of  the  LC.  At 
0  dB,  the  ARLO  detector  is  registering  a  BER  of  10^®'^,  or  0.04  bit  errors  for 
every  100  transmitted.  The  other  detectors’  performance  falls  in  the  range  of 
10^^  ®to0  BER.  While  ARLO  is  outperformed  in  this  scenario,  it  still  provides 
excellent  detection  capabilities. 
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Detector  Performance  in  PB  Interference  Environment 


Figure  3.7:  BER  in  PB  interference  environment. 

3.2.3  Dual  Continuous  Wave  Interferer 

The  2CW  jammer  provides  a  greater  hindrance  to  communication  over  a  channel 
than  does  the  CW  jammer,  as  can  be  seen  from  the  greater  frequency  coverage  of 
the  interferer  in  Figure  3.12,  compared  with  that  of  the  CW  interferer  in  Figure 
3.10.  Therefore,  one  would  expect  that  BER  performance  would  decrease  as 
compared  to  the  single  CW.  For  all  five  detectors,  this  is  the  case,  with  the  four 
non-ARLO  detectors  registering  a  0  dB  BER  in  the  10^^  to  10^5  range,  meaning 
10  —  32  bit  errors  for  every  100  bits  transmitted.  The  ARLO  detector’s  0  dB 
performance  registered  at  slightly  better  than  10^^  BER,  or  less  than  1  error 
out  of  every  100  bits.  ARLO  was  clearly  the  best  performer  in  this  scenario. 

3.2.4  Mixed  Interferer 

The  MX  interference  scenario  is  potentially  the  most  hostile  environment  in  the 
simulation.  While  the  ISR  remains  constant  in  each  environment,  the  frequency 
band  coverage  of  the  MX  interferer  is  greater  than  that  of  the  other  3  jam¬ 
mers.  ARLO  again  clearly  outperforms  the  other  detectors,  registering  a  BER 
of  10^^  ®,  or  0.25  bit  errors  for  every  100  bits  transmitted,  at  the  0  dB  SNR 
point.  The  nearest  other  detector  experiences  10  errors  for  every  100  bits. 
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Power  Spectra  for  PB  Interference  Environment 


Figure  3.8:  Power  spectra  for  PB  interference  environment. 

3.2.5  No  Channel  Interference 

The  purpose  of  running  a  simulation  with  no  interferer  present  is  to  document 
the  effectiveness  of  the  ARLO  detector  in  a  very  stable  environment.  This 
represents  using  the  ARLO  detector  in  a  benign  scenario.  Theory  predicts  that 
the  linear  correlator  will  perform  best,  as  it  is  the  optimum  detector  in  additive 
white  Gaussian  noise.  The  BER  curve  in  Figure  3.15  agrees  with  theory,  with  the 
LC  detector  performing  the  best.  ARLO  performs  only  slightly  less  effectively, 
achieving  a  zero  BER  at  all  SNR  greater  than  (  —  1)  dB.  These  results  indicate 
that  the  ARLO  detector  can  be  effectively  used  even  when  jamming  is  not 
present. 

3.2.6  Performance  Summary 

Figure  3.16,  a  summary  of  section  3.2,  shows  that  the  ARLO  detector  clearly 
outperforms  the  other  three  detectors  in  three  of  the  four  interference  environ¬ 
ments  tested.  In  the  fourth  environment,  consisting  of  the  CW  interferer,  the 
best  detector  achieves  zero  bit  errors,  while  ARLO  performs  admirably  with 
only  0.04  bit  errors  for  every  100  bits  transmitted.  Additionally,  the  ARLO 
detector  never  exceeds  1  bit  error  for  every  100  transmitted  bits,  while  the  best 
of  the  other  detectors  makes  10  bit  errors  75%  of  the  time.  These  results  clearly 
indicate  the  performance  enhancement  that  the  ARLO  detector  can  provide 
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Detector  Performance  in  CW  Interference  Environment 


Figure  3.9:  BER  in  CW  interference  environment. 


Power  Spectra  for  CW  Interference  Environment 


Figure  3.10:  Power  spectra  for  CW  interference  environment. 
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Detector  Performance  in  2CW  Interference  Environment 


Figure  3.11:  BER  in  2CW  interference  environment. 


Pow/er  Spectra  for  2CW  Interference  Environment 


Figure  3.12:  Power  spectra  for  2CW  interference  environment. 
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Detector  Performance  in  MX  Interference  Environment 


Figure  3.13:  BER  in  MX  interference  environment. 


over  other  standard  detectors. 


3.3  Performance  Effects  of  Alternative  Parame¬ 
ter  Values 

A  number  of  parameters  exist  within  the  ARLO  detector  whose  values  are  a 
major  factor  in  its  BER  performance.  Two  such  parameters  were  varied,  to 
determine  their  specific  effect  on  the  ARLO  detector.  Bh,  the  pdf  estimation 
parameter,  defines  the  the  number  of  symbols  used  in  calculating  the  histogram 
used  in  estimating  the  pdf  of  the  received  signal.  P,  the  AR-model  parameter, 
defines  the  number  of  AR  coefficients  used  in  the  AR-model.  In  addition,  the 
ISR,  a  variable  external  to  the  communication  system,  was  varied  for  complete¬ 
ness,  as  the  power  of  an  interference  generator  is  an  unknown  and  can  be  set  to 
any  practical  level. 

3.3.1  Bh:  pdf  Estimation  Parameter 

The  first  variable  tested  was  Bh,  the  number  of  symbols  used  to  calculate  the 
histogram  for  pdf  estimation.  As  will  be  discussed  in  chapter  4,  Bh  has  a  major 
effect  on  the  simulation  runtime.  Therefore,  a  practical  limit  of  Bh  =  128  was 
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Power  Spectra  for  MX  Interference  Environment 


Figure  3.14:  Power  spectra  for  MX  interference  environment. 


Detector  Performance  with  No  Interferer  Present 


Figure  3.15:  BER  with  no  interferer  present. 
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Figure  3.16:  Number  of  bit  errors  at  the  receiver  for  every  100  bits  transmitted, 
at  0  dB  SNR. 


chosen.  A  lower  value  was  also  decided  upon,  Bh  =  64.  Figure  3.17  shows 
the  BER  curve  for  a  mixed  interferer  with  Bh=128.  Comparing  this  with  the 
BER  curve  in  Eigure  3.18,  a  mixed  interference  environment  with  Bh=64,  it  is 
clear  that  a  lower  value  of  Bh  results  in  a  less-smooth  plot.  The  BER  varies 
more  dramatically  with  a  lower  Bh,  making  its  results  less  predictable.  While 
the  results  are  still  good,  the  BER  using  Bh=64  never  improves  beyond  10^®, 
implying  that  a  lower  Bh  imposes  a  performance  limit  on  the  ARLO  detector. 
This  agrees  with  theory,  as  pdf  estimation  improves  when  calculated  using  a 
greater  number  of  symbols.  The  better  the  noise  pdf  is  estimated,  the  greater 
effect  the  ARLO  detector  can  have  in  overcoming  the  noise.  These  results 
are  corroborated  in  Eigures  3.19  and  3.20,  where  a  variable  Bh  is  seen  in  a  CW 
interference  environment.  The  same  unpredictability  is  seen,  with  a  higher  BER 
lower  bound. 
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Detector  Performance  in  MX  Interference  Environment 


Figure  3.17:  BER  curve  in  Mixed  interference  environment  with  Bh  =  128. 


Detector  Performance  in  MX  Interference  Environment 


Figure  3.18:  BER  curve  in  Mixed  interference  environment  with  Bh  =  64. 
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Detector  Perlormance  in  CW  Interference  Environment 
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Figure  3.19:  BER  curve  in  CW  interference  environment  with  Bh  =  128. 


Detector  Performance  in  CW  Interference  Environment 
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Figure  3.20:  BER  curve  in  CW  interference  environment  with  Bh  =  64. 
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3.3.2  P:  AR-model  Parameter 

The  second  variable  tested  was  P,  the  number  of  AR  coefficients  used  in  the  AR- 
model.  P  also  has  an  effect  on  runtime,  as  discussed  in  chapter  4,  so  the  values 
chosen  for  testing  include  P  =  [8,  16,  32] .  Figures  3.21,  3.22,  and  3.23  show 
BER  curves  with  2CW  interferers  and  a  variable  P.  It  is  interesting  to  note  that 
as  P  decreases,  the  performance  of  ARLO  increases.  With  P=32,  the  LOW  and 
LOW  actually  outperform  ARLO,  whereas  with  P=8,  the  ARLO  probability  of 
bit  error  drops  below  10^®.  The  BER  curves,  with  the  CW  interferer,  in  Eigures 
3.24,  3.25,  and  3.26  concur.  The  ARLO  BER  with  P=8  actually  drops  below 
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Detector  Performance  in  2CW  Interference  Environment 


Signal-to-Thermal  Noise  Ratio  (dB) 


Eigure  3.21:  BER  curve  for  2CW  interference  environment  with  P=32. 


3.3.3  ISR:  Interference  Parameter 

Varying  the  ISR  increases  the  practical  significance  of  the  simulations,  as  there 
is  no  guarantee  how  much  interference  power  a  communication  system  may  face 
in  the  channel.  Since  varying  the  ISR  has  no  effect  on  the  simulation  runtime, 
any  practical  values  may  be  chosen.  Tests  were  run  using  ISR  =  [20,  30,  40] . 
Eigures  3.27,  3.28,  and  3.29  show  the  BER  performance  in  a  PB  interference 
environment  with  varying  ISR.  As  expected,  the  performance  of  the  detectors 
falls  when  more  interference  is  present.  ARLO  continues  to  remain  the  best¬ 
performing  detector.  In  Eigures  3.30,  3.31,  and  3.32,  where  the  CW  interferer  is 
examined,  detector  performance  also  decreases  as  ISR  increases.  One  thing  to 
note  is  that  with  ISR=40  dBW,  the  ARLO  detector  begins  outperforming  two 
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Detector  Performance  in  2CW  Interference  Environment 
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Figure  3.22:  BER  curve  for  2CW  interference  environment  with  P=16. 


Detector  Performance  in  2CW  Interference  Environment 
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Figure  3.23:  BER  curve  for  2CW  interference  environment  with  P=8. 
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Detector  Performance  in  CW  Interference  Environment 
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Figure  3.24:  BER  curve  for  CW  interference  environment  with  P=32. 


Detector  Performance  in  CW  Interference  Environment 
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Figure  3.25:  BER  curve  for  CW  interference  environment  with  P=16. 
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Detector  Performance  in  CW  Interference  Environment 
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Figure  3.26:  BER  curve  for  CW  interference  environment  with  P=8. 

other  detectors.  This  shows  the  impressive  robust  nature  of  the  ARLO  detector, 
even  under  heavy  interference. 
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Figure  3.27:  BER  curve  for  PB  interference  environment  with  ISR=20  dbW. 


27 


Detector  Performance  in  PB  Interference  Environment 
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Figure  3.28:  BER  curve  for  PB  interference  environment  with  ISR=30  dbW. 
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Figure  3.29:  BER  curve  for  PB  interference  environment  with  ISR=40  dbW. 
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Detector  Performance  in  CW  Interference  Environment 
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Figure  3.30:  BER  curve  for  CW  interference  environment  with  ISR=20  dbW. 
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Figure  3.31:  BER  curve  for  CW  interference  environment  with  ISR=30  dbW. 
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Detector  Performance  in  CW  Interference  Environment 


Figure  3.32:  BER  curve  for  CW  interference  environment  with  ISR=40  dbW. 

3.4  Autoregressive  Modeling 

Two  methods  were  examined  in  order  to  create  an  autoregressive  model  of  the 
received  signal  vector.  The  Modified  Covariance  Algorithm  (MCA)  was  initially 
chosen  in  this  capacity  because  of  its  simplicity.  The  MCA  generates  a 
order  all-pole  filter  model  to  represent  an  input  signal  by  minimizing  the  sum 
of  the  forward  and  backward  prediction  errors  to  create  a  lattice  model  of  a 
random  process.  The  forward  prediction  error  is  defined  as, 

p 

e'^{n)  =  x{n)  +  ^  a{k)x{n  —  k),  (3.6) 

k=l 


and  the  backward  prediction  error  is, 


p 

e^{n)  =  x{n  —  P)  +  ^  a*{k)x{n  —  P  -\-k),  (3.7) 

k=l 

where  a*  indicates  conjugation,  and  P  is  the  order  of  the  all-pole  filter.  The 
sum  to  be  minimized  is  then, 

JV 

£^=  ^  [|e+(n)|"  +  |e-(n)|"]  .  (3.8) 

n=P 
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To  find  the  filter  coefficients,  a{k),  necessary  to  minimize  ,  (3.8)  is  manipu¬ 
lated,  with  the  result  being  the  normal  equations, 

Ef=i  [fx  ih  k)  +  r^{P  -  k,P  -  /)]  a{k)  =  -  [r^  {1, 0)  -h  (P,  P  -  /)]  qx 

/  =  1,...,P  ’  ^ 

where  rx{l,  k)  is  the  auto-correlation  of  the  process  x{n), 

N 

rx{l,  k)  =  ^  x{n  —  k)x*{n  —  1).  (3.10) 

n=P 

Using  (3.10)  to  obtain  the  auto-covariance  sequences  for  the  random  process 
within  (3.9)  allows  the  normal  equations  to  be  solved  to  obtain  the  filter  coeffi¬ 
cients,  a  {k) . 

The  Burg  method  was  also  examined  due  to  its  performance  characteristics. 
However,  tests  revealed  no  significant  performance  enhancement.  Therefore,  the 
MCA  was  chosen  because  of  its  simplicity. 
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Chapter  4 

Simulation  Conversion 


The  conversion  of  the  simulation  into  a  high-level  language  is  a  major  component 
of  this  research  project,  for  two  reasons.  First,  it  is  expected  that  porting 
the  simulation  into  a  compiled  language  will  increase  the  execution  speed  over 
the  interpreted  code  used  by  MATLAB®  .  Second,  with  the  ultimate  goal 
of  implementing  the  ARLO  detector  on  a  DSP,  many  utilities  exist  that  can 
easily  translate  code  written  in  certain  high-level  languages  into  DSP  machine 
language.  Therefore,  the  porting  of  the  code  is  the  first  step  towards  a  hardware 
implementation  of  the  ARLO  detector. 

This  chapter  will  thoroughly  discuss  the  conversion  process.  The  first  sec¬ 
tion  will  detail  the  runtime  of  the  MATLAB®  simulation,  including  the  effect 
that  certain  parameters  have  on  execution  time.  The  second  section  will  dis¬ 
cuss  using  the  MATLAB®  Compiler  to  port  the  simulation  into  C.  The  design 
and  planning  for  the  conversion  is  covered  in  the  third  section,  including  choice 
of  language  and  platform  options.  The  fourth  section  will  discuss  the  functions 
requiring  conversion,  as  well  as  some  major  obstacles  faced  in  the  conversion  pro¬ 
cess.  Finally,  the  fifth  section  discusses  the  results  of  the  conversion,  and  com¬ 
pares  the  output  of  the  translated  code  to  that  of  the  MATLAB®  simulation. 

4.1  Timing  Analysis 

The  runtime  of  the  MATLAB®  simulation,  defined  as  the  elapsed  time  between 
the  initial  start  time  and  the  time  the  simulation  finishes,  is  a  major  concern 
because  it  is  the  sole  determining  factor  of  whether  the  ARLO  detector  can  be 
used  in  a  real-time  application,  or  whether  its  capabilities  must  be  taken  advan¬ 
tage  of  offline.  The  initial  timing  analysis  was  performed  using  old  hardware, 
with  the  fastest  system  being  a  Pentium  II  300  MHz,  with  256  MB  of  RAM. 
However,  more  current  hardware  has  been  obtained,  and  a  new  timing  analysis 
performed  that  agrees  in  relative  terms  with  the  initial  analysis,  but  with  faster 
overall  execution  speed.  The  new  system  is  a  dual-processor  Pentium  III  XEON 
600MHz,  with  1  GB  of  RAM.  The  results  detailed  below  were  measured  from 
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this  new  system,  using  the  default  parameters  as  defined  in  Table  3.1,  unless 
otherwise  noted. 

All  of  the  timing  results  were  very  consistent,  so  two  simulations  were  chosen 
as  a  representative  sample,  the  CW  and  PB  cases.  The  overall  runtime  for  the 
two  simulations  are  shown  in  Table  4.1,  with  a  difference  in  overall  runtime  of 
1  minute  between  the  two  simulations.  Many  events  comprise  the  fifteen  hour 
runtime  of  the  simulation.  These  include  creating  a  data  stream  to  transmit, 
generating  thermal  background  noise,  generating  the  interferer,  corrupting  the 
signal  in  the  channel,  and  performing  detection  using  all  five  detectors.  This 
being  the  case,  the  fifteen  hour  runtime  is  not,  by  itself,  representative  of  the 
time  needed  for  ARLO  to  perform. 


Simulation  Type 

Runtime  (hrs) 

CW 

15.1119 

PB 

15.1338 

Table  4.1:  Total  runtime,  in  hours,  of  two  representative  simulations. 


Simulation  Performance:  SNR  Runtime 


Figure  4.1:  Runtime  of  simulation  with  CW  interferer,  over  SNR. 

To  obtain  a  greater  understanding  of  the  runtime.  Figures  4.1  and  4.2  show 
the  runtime  of  each  simulation,  broken  down  by  SNR  iteration.  The  variability 
of  the  plots  is  misleading,  however,  as  the  scale  of  the  ordinate  is  expanded.  The 
runtime  of  each  SNR  ratio,  as  shown  in  Figure  4.1,  varies  from  0.7188  hours  to 
0.7202  hours,  a  negligible  difference  of  5  seconds.  For  analysis,  the  average  time 
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Simulation  Performance:  SNR  Runtime 


Figure  4.2:  Runtime  of  simulation  with  PB  interferer,  over  SNR. 


per  SNR  iteration,  for  both  simulations,  is  0.72  hours,  or  43.2  minutes.  Each 
SNR  iteration  is  composed  of  10  Monte  Carlo  (MC)  iterations.  So  the  average 
MC  iteration  runtime  is, 

43.2mm  4.32mm 

- =  - .  (4.1) 

10  iteations  iteration 

More  detailed  study  needs  to  be  performed  to  get  exact  values  (see  Chapter  5) , 
but  the  ARLO  detector,  on  a  single  data  stream  of  16,384  bits,  requires  approx¬ 
imately  four  minutes  to  completely  detect  the  received  data  stream.  While  fast, 
ARLO  is  clearly  not  real-time  capable  in  its  current  incarnation. 

4.1.1  Simulation  Runtime  Effects  of  Alternate  Parameter 
Values 

As  discussed  in  section  3.3,  where  changing  certain  communication  system  pa¬ 
rameter  values  results  in  different  BER  performance  results,  changing  these 
same  parameters  also  often  has  an  effect  on  the  runtime  of  the  simulation.  Table 
4.2  details  the  overall  simulation  runtime  when  changing  the  listed  parameters. 
Changing  the  two  internal  parameters  has  a  noticeable  effect  on  simulation  run¬ 
time.  Halving  Bh  results  in  halving  the  simulation  runtime,  a  clear  1:1  direct 
relationship.  Assuming  all  else  remains  constant,  the  ARLO  detector  takes  ap¬ 
proximately  two  minutes  to  decode  a  16,384  bit  data  stream  with  Bh=64.  This 
is  a  great  runtime  performance  improvement.  But  the  drawback  is  that  the 
BER  performance  degrades,  as  shown  in  section  3.3.1.  Changing  the  value  of 
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Parameters 

CW  Case  (hrs) 

Default  (Bh=128,  P=16,  ISR=30  dB) 

15.1119 

Bh=64 

7.53 

P=8 

13.6932 

II 

CO 

to 

20.2771 

ISR=20  dB 

15.1015 

ISR=40  dB 

15.1235 

Table  4.2:  Simulation  runtime  using  alternate  parameter  values. 


P  also  has  an  effect  on  performance.  Reducing  P  to  8  results  in  a  9.4%  run¬ 
time  performance  increase,  while  increasing  P  to  32  results  in  a  34%  runtime 
performance  loss.  Changing  the  ISR,  a  variable  external  to  the  communication 
system,  has  no  effect  on  runtime,  as  expected. 


4.2  The  MATLAB®  Compiler 

As  discussed  in  the  introduction  to  this  chapter  on  page  32,  it  is  well-known  that 
MATLAB®  code,  while  convenient,  is  inherently  slow  because  the  MATLAB® 
program  is  an  interpreter,  which  directly  executes  high-level  language  code  with¬ 
out  first  compiling  it  into  machine  language.  The  runtime  of  the  simulation  can 
thus  be  greatly  reduced  if  the  mathematically  intensive  functions  of  the  simula¬ 
tion  execute  in  C  or  C++,  languages  whose  code  is  compiled  before  execution. 

The  ideal  solution  is  to  recode  the  simulation  by  hand  into  C  or  C++, 
thereby  optimizing  it  as  much  as  possible.  However,  in  order  to  get  results 
in  the  short  term,  an  automatic,  machine-generated  translation  is  a  good  al¬ 
ternative.  The  MATLAB®  Compiler,  in  conjunction  with  a  standard  compiler 
such  as  Microsoft  Visual  C++,  translates  most  MATLAB®  code  into  C.  It  then 
compiles  the  C  code  into  a  MEX-file,  a  binary  format  directly  accessible  from 
an  M-file  (the  standard  file  used  for  MATLAB®  code),  i.e.  the  M-file  syntax 
to  call  a  function  located  in  another  file  is  the  same  whether  the  target  file  is  a 
MEX-file  or  an  M-file. 

The  major  benefit  of  using  MEX-files  is  its  handling  of  loops.  C  and  C++ 
are  extremely  efficient  in  executing  any  type  of  loop,  whereas  MATLAB®  is 
notoriously  poor  [6,  pp.  4],  simply  due  to  its  status  as  an  interpreter.  This,  in 
itself,  is  not  enough  to  give  a  dramatic  improvement  in  simulation  runtime  be¬ 
cause  the  current  code  uses  few  large  loops.  However,  the  MATLAB®  Compiler 
provides  numerous  compilation  options  that  can  further  enhance  the  runtime. 
They  allow  the  compiler  to  generate  simpler  data  types,  resulting  in  simpler, 
and  thus  faster,  code. 

Two  of  the  options  available  are  particularly  applicable,  and  provide  a  huge 
performance  increase  when  used  in  combination.  The  first  is  the  assertion  that 
all  floating-point  variables  contain  real  numbers.  This  allows  the  compiler  to 
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remove  all  code  whose  purpose  is  to  handle  complex  numbers,  thereby  reducing 
the  complexity  of  the  code  and,  thus,  the  runtime.  The  second  option  generates 
code  that  restricts  matrices  from  growing  beyond  their  initial  size  and,  conse¬ 
quently,  never  checks  a  matrices’  bounds.  Using  these  two  assertions  to  create 
MEX-files,  the  simulation  runtime  decreases  by  roughly  a  factor  of  two. 

One  problem  arose  while  running  simulations  using  the  MEX-files.  Depend¬ 
ing  on  the  specific  parameters  of  the  simulation,  the  pdf  estimation  sometimes 
resulted  in  complex  AR  coefficients.  This  caused  the  simulation  to  stop  exe¬ 
cuting,  with  an  appropriate  error  message.  Recompiling  the  MEX-files  without 
the  real  number  assertion  will  remove  this  error.  However,  that  assertion  is  the 
primary  factor  in  the  speed  increase  of  the  simulation.  Thus,  removing  this 
assertion  will  likely  negate  much  of  the  performance  improvement  [6,  pp.  15]. 

As  promising  as  using  the  MATLAB®  Compiler  appears,  we  receive  no 
performance  benefit  from  it  if  the  speed-increasing  assertions  cannot  be  used. 
Therefore,  to  convert  the  simulation  to  a  high-level  language,  porting  it  by  hand 
is  necessary. 


4.3  Software  Design 

As  in  all  programming  projects,  detailed  planning  must  precede  the  actual  cod¬ 
ing  process.  The  first  important  issue  is  to  decide  which  programming  language 
will  be  optimal  for  this  project.  Both  C  and  C++  were  considered  due  to  their 
common  usage  as  DSP  languages.  Since  C  supports  structured  programming, 
and  this  method  provides  a  good  representation  of  the  simulation,  our  program 
will  not  benefit  from  the  object-oriented  extensions  that  C++  provides.  There¬ 
fore,  C  was  chosen  as  the  programming  language. 

The  next  important  issue  is  the  platform  on  which  the  simulation  will  be  de¬ 
veloped.  One  option  is  Linux,  which  provides  solid  stability,  as  well  as  improved 
performance  and  memory  management  over  Microsoft  Windows®  .  Windows, 
however,  is  more  accessible  to  the  people  who  will  be  performing  the  code  trans¬ 
lation,  and  it  provides  an  enhanced  development  and  debugging  environment  in 
Microsoft  Visual  C++  6.0®  .  Therefore,  Windows  was  chosen  as  the  develop¬ 
ment  platform.  However,  in  order  not  to  restrict  future  development,  we  decided 
to  write  the  code  using  only  standard  ANSI  C  features,  which  should  make  it 
highly  portable  to  Linux,  if  desired. 

The  last  issue  is  to  decide  which  data  type  would  best  represent  the  signal 
samples  used  in  the  simulation.  MIX  Software’s  C/Math  Toolchest©,  which 
provides  vector  and  matrix  math  libraries,  can  be  compiled  to  use  either  float¬ 
ing  point  precision  or  double  precision.  Double  precision  is  clearly  superior  in 
its  ability  to  provide  a  very  accurate  representation  of  the  signals.  But  since 
MATLAB®  uses  double  precision,  and  numerous  problems  have  arisen  regard¬ 
ing  memory  usage  in  past  simulations,  floating  point  precision  was  chosen  to 
determine  if  the  added  precision  is  necessary  or  if  it  poses  an  inefficient  use  of 
computing  resources. 


36 


4.3.1  Simulation  Block  Diagram 

To  convert  the  simulation  into  C,  it  is  imperative  to  understand  the  commu¬ 
nication  system  being  simulated.  An  excellent  way  of  accomplishing  this  is  to 
study  a  visual  representation  of  the  communication  system.  The  block  diagram, 
shown  in  Figure  4.3,  was  created  for  this  purpose. 


Nonlinearity 


(a) 


Estimate  h' 


AR 

Filter 


a  parameters  - 


pdf  Estimate 


►  est[h'] 


(b) 


Figure  4.3:  (a)  ARLO  system  block  diagram,  (b)  Expanded  view  of  Estimate 
h’  block. 

Figure  4.3a  provides  the  overall  signal  flow  of  the  receiver.  Figure  4.3b 
illustrates  the  detail  of  the  nonlinearity  for  estimation  of  the  function  h’  (to 
be  described  later).  In  Figure  4.3a,  p  is  the  received  signal  vector  of  length 
N,  which  is  composed  of  both  the  signal  sent  by  the  transmitter  and  the  noise 
present  in  the  channel.  The  receiver  must  first  calculate  the  LO  nonlinearity, 
g  {p) ,  which  is  individually  calculated  for  each  incoming  symbol  in  p,  each  of 
which  is  composed  of  samples  per  symbol.  The  next  step  is  to  de-spread  the 
“chipped”  data  by  using  the  pseudo-noise  (PN)  sequence.  Finally,  an  estimate 
of  the  transmitted  signal,  est[s],  is  computed  by  summing  the  de-spread  samples 
of  the  symbols  and  thresholding. 

The  receiver  output,  est[s],  is  highly  dependent  on  the  nonlinearity  compu¬ 
tation.  This  computation  is  a  multi-step  process.  The  first  step  in  calculating 
the  LO  nonlinearity  is  the  computation  of  the  P*^-order  AR  model  coefficients. 


37 


{di} ,  i  =  0,1,2, P.  To  date,  the  modified  covariance  algorithm  (MCA)  is 
used  for  this  purpose.  The  second  step  is  to  estimate  the  pdf  of  the  noise.  This 
is  currently  accomplished  via  a  histogram  approximation  method.  This  process 
is  shown  in  detail  in  Figure  4.3b.  The  received  signal,  p,  passes  through  the 
AR  filter,  with  coefficients  {di} ,  as  shown  in  Figure  4.3a,  to  create  a  whitened 
signal  vector,  w.  The  histogram  of  w  is  then  computed.  The  result  represents 
the  pdf  of  the  whitened  received  signal 

The  pdf  of  the  noise  component  of  w  can  be  written  as. 


fniv)  =  fn{Vl,V2,---,VN) 

=  UZl  fnivil  Vi-1,  ■■■,Vi-p)  ' 

where 


frii  {Vi  I  Vi—1 ,  •  •  •  ,  Vi—p) 


This  can  be  rewritten  as. 


fni  iVi) ,  fori  =  1 
fm  ivi  I  Vi-1,- ■■,Vi), 
for  i  =  2, . . . ,  P 


fniv) 


N  (  ^ 

fw  I  ~  ^jVi-j 

i=l  \  j=0 


(4.2) 


(4.3) 


(4.4) 


where  (w)  is  the  pdf  of  the  white  noise  process  and  oq  =  —1.  The  result  of 
this  block,  as  estimate  of  h',  est[h’],  is  the  derivative  of  the  natural  log  of  the 
pdf  of  w,  that  is, 

est  [h'  (w)]  =  ^  {In  [f^  (w)]}  .  (4.5) 

The  next,  and  final,  step  in  calculating  g  (p)  is  to  use  p  and  h’  in  the  filtering 
process. 


niin(P,iV— i)  /  p 


9i  i'u)  = 

E 

’  for  i  =  1, . . 

..,N,  (4.6) 

V  ^=0  ) 

where 

J  0,  for  i  <  0 

(4.7) 

Vi 

“  1  Pi,  for  i  G  [1,  A^] 

4.4  Code  Translation 

The  simulation  is  composed  of  two  distinct  sections.  The  first  section  imple¬ 
ments  the  direct  sequence  spread  spectrum  generator,  the  communication  chan¬ 
nel,  including  a  thermal  noise  generator  and  interference  generator,  and  the 
linear  correlator  receiver.  The  second  part  of  the  simulation  contains  the  non¬ 
linear  LO  receivers.  The  conversion  into  C  was  performed  along  these  lines. 
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4.4.1  Function  List 

Since  the  research  team  decided  to  port  the  simulation,  based  on  the  MATLAB® 
simulation,  into  C,  the  functional  breakdown  required  in  C  will  be  very  similar 
to  that  used  in  MATLAB®  .  The  list  of  functions  in  the  completed  C  code 
include  the  following: 

•  ar_histo.c:  Estimates  the  pdf  of  a  white  noise  process,  without  assuming 
any  underlying  structure  for  the  white  noise  pdf. 

•  ar_theory.c:  Estimates  the  pdf  of  a  white  noise  process,  assuming  that 
the  pdf  of  the  white  noise  has  a  known  underlying  structure,  i.e.  Gaussian, 
Poisson,  or  Laplacian. 

•  contwave.c:  Creates  CW  interferer. 

•  conv.c:  Performs  convolution  of  two  vectors  of  polynomial  coefficients. 

•  decision. c:  Makes  decision  on  received  symbol. 

•  dsss.c:  Generates  direct-sequence  spread-spectrum  signal. 

•  fftshift.c:  Swaps  left  and  right  halves  of  an  input  vector,  placing  the  DC 
(frequency=0)  component  in  the  middle  of  the  spectrum. 

•  filter. c:  Filters  a  data  stream  based  on  input  filter  coefficients. 

•  fliplr.c:  Flips  the  columns  of  a  matrix  in  the  left-right  direction  about  a 
vertical  axis. 

•  flipud.c:  Flips  the  rows  of  a  matrix  in  the  up-down  direction  about  a 
horizontal  axis. 

•  hist.c:  Creates  histogram  of  a  data  stream. 

•  modcov.c:  Finds  P*^-order  all-pole  model  for  a  signal  using  the  MCA 
method. 

•  nlp.c:  Calculates  the  LO  nonlinearity,  g  (p) . 

•  read_parms.c:  Inputs  the  simulation  parameters  from  a  configuration  file. 

•  receiver. c:  Performs  synchronous  reception  of  a  transmitted  signal. 

•  reshape. c:  Reshapes  a  matrix. 

•  simulation. c:  Main  program  simulation  file. 

•  th_inter.c:  Creates  a  theoretical  PB  interferer. 

•  toeplitz.c:  Creates  a  Toeplitz  matrix  from  two  input  vectors. 

The  list  of  header  files  include: 
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•  mathlib.h:  Header  file  for  C/Math  Toolchest®  . 

•  function_list.h:  Function  prototypes. 

•  sim_parameters.h:  Simulation  global  variable  declarations. 

•  sim_parameters_extern.h:  Simulation  global  variable  declarations,  exter¬ 
nalized. 

The  complete  source  code  for  the  simulation  can  be  found  in  Appendix  A. 

4.4.2  Plotting  Scripts 

After  debugging  was  complete,  the  research  team  required  a  method  to  plot 
the  results  of  the  simulation,  in  order  to  compare  its  output  with  that  of 
MATLAB®  .  With  the  final  goal  of  running  the  ARLO  detector  on  a  DSP,  it  is 
unnecessary  to  write  a  routine  in  C  which  will  handle  BER  plotting,  especially 
since  graphical  coding  in  C  would  add  a  level  of  complexity  to  the  conversion 
attempt  that  is  not  required.  Since  MATLAB®  is  particularly  suited  for  graph¬ 
ical  plotting,  a  script  was  written  in  MATLAB®  that  can  read  in  the  output  of 
a  C  simulation  from  a  file  and  graph  the  BER  and  spectral  plots.  These  plots 
can  then  be  compared  with  those  from  a  MATLAB®  simulation,  for  testing 
purposes.  The  source  code  for  this  script  can  be  seen  in  Appendix  B. 


4.5  Conversion  Results 

The  results  of  the  conversion  were  highly  successful.  Comparisons  of  BER  per¬ 
formance,  using  identical  system  parameters,  of  the  C  and  MATLAB®  simu¬ 
lations  are  shown  in  Eigures  4.4,  4.5,  4.6,  and  4.7.  Erom  these  figures,  it  is  clear 
that  the  simulation  has  been  successfully  ported  to  C,  as  the  BER  curves  are 
nearly  identical. 

The  power  spectra  curve  of  the  PB  scenario,  shown  in  Eigures  4.8  and  4.9, 
are  also  nearly  identical,  therefore  showing  that  the  conversion  was  a  success. 


40 


Detector  Pertormance  in  PB  Interference  Environment 


Signal-to-Thermal  Noise  Ratio  (dB) 


Figure  4.4:  C  Simulation:  BER  curve  in  PB  interference  environment. 


Detector  Performance  in  PB  Interference  Environment 


Signal-to-Thermal  Noise  Ratio  (dB) 


Figure  4.5:  MATLAB®  Simulation:  BER  curve  in  PB  interference  environment. 
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Detector  Performance  in  2CW  Interference  Environment 


Figure  4.6:  C  Simulation:  BER  curve  in  2CW  interference  environment. 


Detector  Performance  in  2CW  Interference  Environment 


Figure  4.7:  MATLAB®  Simulation:  BER  curve  in  2CW  interference  environ¬ 
ment. 
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Power  Spectra  for  PB  Interference  Environment 


Figure  4.8:  C  Simulation:  Power  spectra  of  PB  interference  environment. 


Power  Spectra  for  PB  Interference  Environment 


Figure  4.9:  MATLAB®  Simulation:  Power  spectra  of  PB  interference  environ¬ 
ment. 
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Chapter  5 

Summary  and  Future 
Research 


In  summary,  the  work  performed  by  the  research  team  was  highly  successful. 
The  performance  tests  expanded  on  previous  research  to  prove  the  efficacy  of 
the  ARLO  detector  in  intense  interference  environments,  as  well  as  its  supe¬ 
rior  BER  performance  in  most  every  scenario  over  the  other  popular  detector 
algorithms.  Simulation  results  also  showed  that,  in  a  non-hostile  channel  with 
zero  interference,  the  performance  differential  between  ARLO  and  the  optimum 
linear  correlator  were  negligible. 

The  conversion  of  the  simulation  into  C  was  also  a  major  milestone  for 
the  research  team.  Intelligent  choices  were  decided  upon  with  regard  to  the 
fundamentals  of  the  conversion,  including  its  platform,  chosen  programming 
language,  and  library  selection.  The  conversion  was  completed  and  thoroughly 
tested  to  be  equivalent  to  its  MATLAB®  predecessor. 

Future  work  is  needed  to  bring  the  ARLO  detector  into  service.  The  C 
simulation  performs  at  a  much  slower  rate  than  the  MATLAB®  version.  Since 
the  C  version  is  necessary  for  downloading  the  detector  onto  a  digital  signal 
processor  (DSP),  a  detailed  study  of  the  required  operating  time  of  each  module 
of  the  simulation  is  required  to  examine  the  feasibility  of  the  ARLO  detector  for 
real-time  applications.  To  accomplish  these  goals,  the  proposed  research  tasks 
include  the  following: 

•  Perform  a  detailed  study  and  comparison  of  the  runtime  of  the  C  and 
MATLAB®  codes 

•  Optimize  the  C  code  to  achieve  maximum  detector  runtime  performance 

•  Research  other  matrix  libraries  for  decreased  runtime 

•  Study  the  advantages  and  disadvantages  of  adapting  the  C  code  to  execute 
on  a  Linux-based  system 
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•  Convert  the  C  code  for  use  on  a  Linux-based  system 

•  Continue  research  on  Parzen’s  estimator,  and  compare  its  performance 
with  that  of  the  histogram  method. 


5.1  Future  Research  Details 

5.1.1  Runtime  Study 

The  first  area  of  future  research  is  to  perform  a  detailed  timing  analysis  of 
the  C-coded  simulation  (CCS).  This  includes  gathering  runtime  data  for  each 
individual  module  of  the  CCS,  comparing  it  to  the  timing  results  previously 
obtained  from  the  MATLAB®  -coded  simulation  (MCS),  and  determining  which 
processes  are  consuming  the  most  execution  time.  The  most  time-consuming 
modules  will  be  examined  more  closely  to  ascertain  the  reasons  behind  this 
performance  lag.  Parameter  variation,  and  its  effect  on  runtime,  will  also  be 
examined  through  structured  simulations  and  timing  comparisons. 

5.1.2  Code  Optimization 

The  current  execution  speed  of  the  CCS  is  below  expectations  due  to  the  re¬ 
search  team’s  previous  focus  on  obtaining  total  functionality  at  the  expense  of 
other  factors.  For  effective,  real-time  application,  it  is  imperative  to  increase 
execution  speed.  One  method  to  accomplish  this  is  to  optimize  the  code  itself, 
focusing  primarily  on  modules  exhibiting  large  execution  time  as  measured  in 
the  runtime  study  (see  5.1.1).  In  addition,  the  compilers  and  linkers  used  in 
building  this  application  have  a  great  many  optimization  options  to  exploit  for 
increased  performance  of  the  CCS.  The  simulation  can  also  be  adapted  to  run 
without  including  debugging  information  in  the  compiled  executable,  which  will 
potentially  result  in  a  noticeable  performance  increase. 

5.1.3  Matrix  Libraries 

To  provide  standard  matrix  and  vector  operations  and  manipulations,  the  CCS 
makes  use  of  the  C/Math  Toolchest'^^  code  library.  This  code  library,  while 
inexpensive,  easy  to  use,  and  fairly  complete,  is  not  optimized  for  current  pro¬ 
cessors.  Another  area  of  future  research  is,  therefore,  to  locate  and  study  other 
matrix  libraries,  both  open-source  and  commercial,  and  determine  their  suit¬ 
ability  for  the  CCS.  If  a  library  is  found  that  can  potentially  enhance  detector 
performance,  the  research  team  will  use  its  extensive  programming  experience 
to  adapt  the  CCS  to  use  the  new  library.  The  team  will  then  perform  a  new  de¬ 
tailed  timing  analysis  to  determine  the  exact  performance  benefits  reaped  from 
its  use. 
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5.1.4  Platform  Migration 

The  CCS  was  coded  using  Microsoft  Visual  C++  6.0,  running  under  Microsoft 
Windows  2000.  It  was  written  using  standard  ANSI  C  (no  Microsoft  extensions 
to  C  were  used)  to  retain  maximum  portability.  The  proposed  research  includes 
examining  whether  CCS  performance  under  Linux  using  a  native  C  compiler, 
such  as  the  GNU  Compiler  Collection  (GCC),  may  give  improved  performance 
over  the  Microsoft-centric  implementation.  In  the  event  that  it  does,  the  re¬ 
search  team  will  migrate  the  code  to  Linux  to  harness  the  performance  benefits 
of  this  operating  system. 

5.1.5  Parzen’s  Estimator 

Further  research  is  needed  in  order  to  determine  if  this  new  method  definitely 
outperforms  the  histogram  method.  The  performance  time,  input  data  size 
requirements,  as  well  as  the  diversity  of  each  method  must  be  more  closely 
examined.  Future  research  may  also  include  testing  the  performance  of  the 
Parzen  estimator  when  alternative  RBFs  are  used,  as  well  as  possibly  researching 
other  pdf  estimation  methods. 
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Appendix  A 


Simulation  Source  Code 

A.l  Header  Files 

A. 1.1  Function  list.h 


/*  Function  prototypes  */ 

void  ar_histo  (Real_Vector  stimulus,  unsigned  stimulus_length,  Real_Vector  a, 

unsigned  a_length,  Real_Vector  h,  Real_Vector  bin_width, 
Real_Vector  mini) ; 

/*  bin_width  &  mini  are  vectors  to  enable  pass-by-reference  */ 
void  ar_theory  (int  iidparmsH,  char  iidtype[],  int  ISR,  Real  K,  Real  T, 

Real_Vector  h) ; 

void  contwave  (Real_Vector  interferer) ; 

void  conv  (Real_Vector  vl,  int  lengthl,  Real_Vector  v2,  int  lengths, 

Real_Vector  output) ; 

void  decision  (Real_Vector  pn,  Real_Vector  g,  Real_Vector  decision); 
void  dsss  (Real_Vector  ref,  int  amp,  Real_Vector  coded.data, 

Real_Vector  chips) ; 

void  fftshift  (Complex_Vector  vl,  unsigned  int  lengthl); 
void  filter  (Real_Vector  b,  int  M,  Real_Vector  a,  int  N, 

Real_Vector  data_in,  int  input_length,  Real_Vector  y) ; 
Real.Natrlx  fllplr  (Real.Natrlx  M,  unsigned  m,  unsigned  n) ; 

Real.Matrlx  fllpud  (Real.Natrlx  M,  unsigned  m,  unsigned  n) ; 

void  hist  (Real.Vector  signal,  unsigned  slgnal.length,  unsigned  nblns, 

Real_Vector  freq,  Real_Vector  bincenter) ; 
void  modcov  (Real.Vector  x,  unsigned  x.length,  unsigned  p, 

Real_Vector  ARparm,  Real  error) ; 

Real_Vector  nip  (Real_Vector  r,  unsigned  int  r_length,  Real_Vector  mini, 

Real_Vector  h,  unsigned  int  h_length,  Real_Vector  bw, 
Real.Vector  a,  unsigned  Int  a.length, 
unsigned  Int  g.length) ; 

Int  read.parms  (char  ^filename); 

void  receiver  (Real_Vector  signal,  Real_Vector  noise,  char  *limit, 

Real_Vector  rcvd_tr) ; 

Real.Matrlx  reshape  (Real.Natrlx  M,  unsigned  ml,  unsigned  nl,  unsigned  m2, 

unsigned  n2) ; 

void  th_inter  (Real_Vector  interferer) ; 

Real_Matrix  toeplitz  (Real_Vector  coll,  Real_Vector  rowl,  unsigned  int  ncl, 

unsigned  Int  nrl); 


A.  1.2  Sim  parameters. h 


/*  Simulation  Parameters 


*  By;  Michael  Banys 

* 

*  January  26,  2000 
*/ 
int 


/*  Book  keeping  parameters  */ 
it,  /*  Number  of  Monte  Carlo  iterations  */ 

Bh,  /*  Total  number  of  symbols  to  calculate  histogram  */ 

K,  /*  Data  points  to  calculate  */ 

T,  /*  Range  of  support  for  pdf  */ 

z,  /*  Multiplier  for  determining  number  of  transmitted  symbols  */ 
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/*  Spread  spectrum  parameters  */ 

Nc,  /*  Samples  per  chip  */ 

PG,  /*  Processing  gain  (chips  per  bit)  */ 

Nb,  /*  Samples  per  symbol  (assume  symbol==bit)  */ 

B,  /*  Symbols  transmitted/received  and  sent  to  BER  tester  */ 
Ntot,  /*  Total  number  of  samples  in  run  */ 

/*  AR  Model  parameters  */ 

Bar,  /*  Total  number  of  symbols  to  estimate  AR  parms  */ 

Bg,  /*  Total  number  of  symbols  to  calculate  nip  */ 

/*  Signal  parameters  ft  Thermal  Noise  */ 
sigmat,  /*  Standard  deviation  of  thermal  noise  */ 
slgmat2,  /*  Variance  of  thermal  noise  */ 

Eb_min,  /*  Minimum  value  of  thermal  SNR  in  dB  */ 

Eb_max,  /*  Maximum  value  of  thermal  SNR  In  dB  */ 

/*  Correlated  interferer  parameters  */ 

ISR,  /*  Interferer-to-Signal  ratio  (dB)  */ 

/*  Continuous  Wave  Interference  */ 

I,  /*  Number  of  CW  interferers:  1  or  2  */ 

P,  /*  Number  of  AR  filter  coefficients  */ 


/*  Wlde-band  filtered  whlte-nolse  jammers  */ 
L,  /*  ...  for  a  filter  of  order  2''L  */ 
iidparms[2],  /*  Parameter  vector  for  iidtype  */ 


/*  Mixed  Jammer  */ 

ISRcw,  /*  Divide  up  interferer  power  */ 
ISRpb; 


/*  Continuous  Wave  Jammer  */ 
freq[2] ,  /*  CW  Interferer  frequency(s)  */ 


polynomial [2] ; 
char 


/*  Wlde-band  filtered  whlte-nolse  jammers  */ 


/*  AR  Model  Parameters  */ 


/*  The  pole  of  the  low-pass  filter  */ 


/*  Coefficient  estimation  method 
’mcov’  =  Modified  Covariance  */ 


limit [5] , 


/*  Signal  Parameters  and  Thermal  Noise  */ 


/*  Limiter  on  noise: 

’hard’  -  noise  blanker 
’soft’  -  soft  limiter  */ 


/*  Correlated  interferer  parameters  */ 

Itype[3] , 


/*  Interference  type: 

’cw’  -  Continuous  Wave 

’pb’  -  Partial  Band 

’mx’  -  Mixed  Partial  Band  and  CW  */ 


/*  Wideband  filtered  whlte-nolse  jammer  */ 

iidtype [3];  /*  The  type  of  lid  noise  sequence  */ 

Real_Vector  Eb_slgma,  poly,  a_t; 


A. 1.3  Sim_parameters_extern.h 

/*  Simulation  Parameters  (extern) 


*  By:  Michael  Banys 

* 

♦  January  26,  2000 
*/ 

extern  Int 

/*  Book  keeping  parameters  */ 
it,  /*  Number  of  Monte  Carlo  iterations  */ 

Bh,  /*  Total  number  of  symbols  to  calculate  histogram  */ 

K,  /*  Data  points  to  calculate  */ 

T,  /*  Range  of  support  for  pdf  */ 

z,  /*  Multiplier  for  determining  number  of  transmitted  symbols  */ 

/*  Spread  spectrum  parameters  */ 

Nc,  /*  Samples  per  chip  */ 

PG,  /*  Processing  gain  (chips  per  bit)  */ 

Nb,  /*  Samples  per  symbol  (assume  symbol==bit)  */ 

B,  /*  Symbols  transmitted/received  and  sent  to  BER  tester  */ 

Ntot,  /*  Total  number  of  samples  in  run  */ 

/*  AR  Model  parameters  */ 

Bar,  /*  Total  number  of  symbols  to  estimate  AR  parms  */ 

Bg,  /*  Total  number  of  symbols  to  calculate  nip  */ 

/*  Signal  parameters  ft  Thermal  Noise  */ 
sigmat,  /*  Standard  deviation  of  thermal  noise  */ 
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slgmat2,  /*  Variance  of  thermal  noise  */ 

Eb_mln,  /*  Minimum  value  of  thermal  SNR  In  dB  */ 
Eb_max,  /*  Maximum  value  of  thermal  SNR  In  dB  */ 

/*  Correlated  Interferer  parameters  */ 
ISR,  /*  Interferer-to-Slgnal  ratio  (dB)  */ 

/*  Continuous  Wave  Interference  */ 

1,  /*  Number  of  CW  Interferers:  1  or  2  ♦/ 

P,  /*  Number  of  AR  filter  coefficients  */ 


/*  Wlde-band  filtered  whlte-nolse  jammers  */ 
L,  /*  ...  for  a  filter  of  order  2''L  */ 
lldparms[2],  /*  Parameter  vector  for  lldtype  */ 

/*  Mixed  Jammer  */ 

ISRcw,  /*  Divide  up  Interferer  power  */ 

ISRpb; 

extern  float 

/*  Continuous  Wave  Jammer  */ 
freq[2] ,  /*  CW  Interferer  frequency(s)  */ 


/* 

polynomial [2] ; 
extern  char 


Wlde-band  filtered  whlte-nolse  jammers  */ 

/*  The  pole  of  the  low-pass  filter  */ 


/*  AR  Model  Parameters  */ 

ar_method[S] , 


/*  Coefficient  estimation  method 
’mcov’  =  Modified  Covariance  */ 


limit [5] , 


/*  Signal  Parameters  and  Thermal  Noise  */ 


/*  Limiter  on  noise: 

’hard’  -  noise  blanker 
’soft’  -  soft  limiter  */ 


/*  Correlated  Interferer  parameters  */ 

Itype[3] , 


/*  Interference  type: 

’cw’  -  Continuous  Wave 

’pb’  -  Partial  Band 

’mx’  -  Mixed  Partial  Band  and  CW  */ 


/* 

lldtype [3] ; 
extern  Real_Vector 


Wideband  filtered  whlte-nolse  jammer 
Eb_slgma,  poly,  a_t; 


*/ 

/*  The  type  of  lid  noise  sequence  */ 


A. 2  Program  Code  Files 


A. 2.1  Simulation. c 

/*  Program  Simulation  main  file 
* 

* 

*  By:  Michael  Banys 
*/ 

#lnclude  <stdlo.h> 

#lnclude  <strlng.h> 

#lnclude  <stdllb.h> 

#lnclude  <tlme.h> 

#lnclude  <math.h> 

#lnclude  "mathllb . h" 

#lnclude  "slm_parameters .h" 

#lnclude  "functlon_llst .h" 
void 

main  (Int  argc,  char  *argv[]) 

i 

/*  PROGRAM  DECLARATIONS  */ 

Int  1,  j , 
saveData, 
where_runnlng, 
snr, 
pe.avg, 
amp, 

Nsamh, 

hlstonum, 

symcnt; 

Real  sym.amp, 
num_errs_lc, 
num_errs_wc, 
num_errs_lo, 
num_errs_lw, 
num_errs_ar, 
error; 

Real_Vector  temp,  temp2, 
ref, 
chips, 
coded. data. 


/*  Counter  variable  */ 

/*  To  control  saving  data  only  once  */ 

/*  To  determine  where  to  save  data  */ 

/*  SNR  Iteration  */ 

/*  Monte  Carlo  Iteration  */ 

/*  Voltage  amplitude  of  signal  */ 

/*  Number  of  samples  In  each  histogram  =  Bh*Nb  */ 

/*  Process  Bh*Nb  samples  at  a  time  for  a  histogram  (loop  var) 
/*  Counter  for  each  symbol  In  AR  section  */ 

/*  Symbol  amplitude  assuming  chip-matched  filtering  */ 

/*  Number  of  errors  for  linear  correlator  */ 

/*  Number  of  errors  for  Ic  on  whitened  */ 

/*  Number  of  errors  for  lo  */ 

/*  Number  of  errors  for  Iw  */ 

/*  Number  of  errors  for  ar  */ 

/*  Error  In  modcovO  */ 


/*  Temp  storage  */ 


/*  Sampled 


/*  BPSK  modulated  signal  */ 
/*  PN  Sequence  */ 

&  PN-modulated  data  signal  */ 


*/ 
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/*  Interferer  */ 


noise,  /*  Total  noise  environment  */ 

rcvd_tr,  /*  Received  signal  (r=s+n)  */ 

whitened,  /*  Whitened  rcvd_tr  */ 

a_hat,  /*  AR  parameter  estimate  vector  */ 

a_hat_negl,  /*  As  above,  except  a_hat_negl [0] =-l,  length  =  1  ♦/ 

g_lo,  g_ar,  g_lw,  /*  Initialize  LO,  AR,  and  LO  on  Wh.  detectors  */ 

g_temp,  /*  Temp  storage  for  g,  holding  one  symbol’s  worth  */ 

lc_dec,  /*  linear  correlator  decision  vector  */ 

lo_dec,  /*  Locally  Optimium  decision  vector  */ 

lw_dec,  /*  LO  on  whitened  decision  vector  */ 

ar_dec,  /*  ARLO  decision  vector  */ 

wc_dec,  /*  Ic  on  whitened  decision  vector  */ 

Pb_lc_avg,  /*  Average  P[err]  */ 

Pb_wc_avg,  Pb_lo_avg,  Pb_lw_avg,  Pb_ar_avg,  block,  /*  Block  to  prepare  for  histogram  */ 

symbol,  /*  One  symbol  */ 

symbolindx,  /*  Index  to  separate  out  symbols  from  bitstream  */ 

filsym,  /*  whitened  symbol  */ 

h_lo,  h_wh,  /♦  LO  h()  ft  AR  h()  ♦/ 

bw_lo,  bw_wh,  /*  LO  ft  AR  bin_width  */ 

min_lo,  mln_wh;  /*  min  value  of  histogram/pdf  LO  ft  AR  */ 

Real_Matrix  Pb_lc,  /*  P[err]  for  Ic  */ 

Pb_wc,  Pb_lo,  Pb_lw,  Pb_ar; 

FILE  *Pblc,  *Pbwc,  *Pblo,  *Pblw,  *Pbar,  *snr_vals, 

*interferer_out,  *noise_out,  *rcvd_tr_out,  *a_hat_out,  *ref_out; 
long  Itime;  /*  For  randO  usage  */ 

if  (argc  !=  2)  /*  Parameter  file  check  */ 

■{ 

printf  ("You  forgot  to  enter  the  parameter  filenameAn") ; 
exit  (1) ; 

} 

read.parms  (argv[l]);  /*  Read  simulation  parameters  from  file  */ 

/*********  Vector  and  matrix  creation  **********/ 

ref  =  valloc  (NULL,  B) ;  /*  Create  signal  */ 

coded.data  =  valloc  (NULL,  Nb  *  B) ;  /*  Sampled  ft  PN-modulated  data  signal  */ 

chips  =  valloc  (NULL,  Nb  *  B) ;  /*  Chip  sequence  */ 

interferer  =  valloc  (NULL,  Ntot) ;  /*  Create  interferer  */ 


chips  =  valloc  (NULL,  Nb  *  B) ; 
interferer  =  valloc  (NULL,  Ntot) ; 
noise  =  valloc  (NULL,  Ntot) ; 
rcvd_tr  =  valloc  (NULL,  Ntot); 
whitened  =  valloc  (NULL,  Ntot) ; 


/*  whitened  rcvd_tr  */ 


a_hat  =  valloc  (NULL,  P  +  1) ;  /♦  AR  parameter  estimate  vector  */ 


/*  block  to  prepare  for  histogram  */ 
/*  Single  symbol  */ 


/*  whitened  symbol  */ 

/*  allocate  decision  vector  */ 


a_hat_negl  =  valloc  (NULL,  1);  /*  a_hat[o]=-l,  length(a_hat)  =  1  */ 

g_lo  =  valloc  (NULL,  Ntot);  /*  LO  detector  */ 

g_ar  =  valloc  (NULL,  Ntot);  /*  AR  detector  */ 

g_lw  =  valloc  (NULL,  Ntot);  /*  LO  on  Wh.  detector  */ 

g_temp  =  valloc  (NULL,  Nb) ;  /*  To  hold  one  symbol’s  worth  of  information  */ 

Pb_lc  =  mxalloc  (NULL,  it,  Eb_max  -  Eb_min  +  1); 

Pb_wc  =  mxalloc  (NULL,  it,  Eb_max  -  Eb_min  +  1); 

Pb_lo  =  mxalloc  (NULL,  it,  Eb_max  -  Eb_min  +  1); 

Pb_lw  =  mxalloc  (NULL,  it,  Eb_max  -  Eb_min  +  1); 

Pb_ar  =  mxalloc  (NULL,  it,  Eb_max  -  Eb_min  +  1); 

Pb_lc_avg  =  valloc  (NULL,  Eb_max  -  Eb_min  +  1); 

Pb_wc_avg  =  valloc  (NULL,  Eb_max  -  Eb_min  +  1); 

Pb_lo_avg  =  valloc  (NULL,  Eb_max  -  Eb_min  +  1); 

Pb_lw_avg  =  valloc  (NULL,  Eb_max  -  Eb_min  +  1); 

Pb_ar_avg  =  valloc  (NULL,  Eb_max  -  Eb_min  +  1); 

block  =  valloc  (NULL,  Bh  *  Nb) ;  /*  block  to  prepare  for  histogram  */ 

symbol  =  valloc  (NULL,  Nb) ;  /*  Single  symbol  */ 

symbolindx  =  valloc  (NULL,  Nb) ; 

filsym  =  valloc  (NULL,  Nb) ;  /*  whitened  symbol  */ 

wc_dec  =  valloc  (NULL,  B) ;  /*  allocate  decision  vector  */ 

lc_dec  =  valloc  (NULL,  B) ; 

lo_dec  =  valloc  (NULL,  B) ; 

lw_dec  =  valloc  (NULL,  B) ; 

ar_dec  =  valloc  (NULL,  B) ; 

h_lo  =  valloc  (NULL,  K  +  l) ;  /*  "h"  function  */ 

h_wh  =  valloc  (NULL,  K  +  l) ; 

bw_lo  =  valloc  (NULL,  l) ;  /*  bin_width  LO  ft  AR  */ 

bw_wh  =  valloc  (NULL,  1); 

min_lo  =  valloc  (NULL,  1) ;  /*  min  value  of  histogram  LO  ft  AR  */ 

min_wh  =  valloc  (NULL,  1) ; 

/t**********  End  creation  section  ***********/ 

a_hat_negl[0]  =  -1.0;  /*  Assign  value  */ 

saveData  =  1;  /*  Save  data  first  time  through  */ 

where .running  =1;  /*  l=Sundevll,  2=L90  */ 

for  (snr  =  0;  snr  <  Eb.max  -  Eb.min  +  1;  snr++)  /*  Iterate  over  each  snr  */ 

■{ 

for  (pe.avg  =  0;  pe.avg  <  it;  pe_avg++)  /*  Iterate  Monte  Carlo  */ 

/*  For  randO  */ 

Itlme  =  time  (NULL); 

srand  ((unsigned)  Itlme  /  2) ;  /*  compute  seed  for  randO  */ 

for  (1=0;  1  <  B;  i++)  /*  Create  random  signal  */ 

■{ 

♦(ref  +  i)  =  (2  ♦  (Real)  rand  0  /  RAND.MAX)  -  1.0; 
if  (♦(ref  +  i)  <  O) 

♦(ref  +  1)  =  -1.0; 
else  if  (♦(ref  +  i)  >  0) 

♦(ref  +  i)  =  1.0; 


/♦  bin.width  LO  ft  AR  ♦/ 


/♦  min  value  of  histogram  LO  ft  AR  ♦/ 


♦(ref  +  i)  =  0.0; 
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} 

syn.anp  =  sqrt  ((double)  slgniat2  *  Eb_slgnia[snr] ) ;  /*  Symbol  amplitude  assuming  chlp- 


/*  NOTE;  Let  the  signal  have  amplitude  of  ’1’  and  then  divide  the  thermal  noise 
by  the  symbol  amplitude  (below,  In  ’Implement  the  channel’  )  */ 
amp  =  1; 

dsss  (ref,  amp,  coded.data,  chips);  /*  Transmitted  DSSS  signal  */ 

prlntf  ("DSSS  done,  [It ,snr]=[Xl,Xl] \n" ,  pe_avg,  snr) ;  /*  Progress  Indicator  */ 

/t**********  Implement  the  Channel  **************/ 
vlnlt  (Interferer,  Ntot,  0.0);  /*  Initialize  Interferer  to  zero  */ 

If  ((ItypeEo]  ==  ’c’)  &&  dtypeCl]  ==  ’w’)) 
contwave  (Interferer) ; 

else  If  ((ItypeEo]  ==  ’p’)  (ItypeCl]  ==  ’b’)) 
th_lnter  (Interferer) ; 

else  If  ((ltype[0]  ==  ’m’)  (Itype[l]  ==  ’x’)) 

■E 

th_lnter  (Interferer) ; 
contwave  (Interferer) ; 

} 

prlntf  ("Interferer  done,  [It ,snr]=[Xl,Xl]\n" ,  pe_avg,  snr);  /*  Progress  Indicator  */ 

Add  thermal  noise  with  a  Gaussian  pdf  of  N(mu=0,  slgmat''2) 
to  the  correlated  Interferer  (and  additive  channel  Is  assumed  )  */ 
for  (1=0;  1  <  Ntot;  1++) 

* (noise  +  1)  =  * (Interferer  +  1)  + 

(((double)  slgmat  /  sym_amp)  *  normal  (1,  0)); 

/t**********  Implement  the  observed  signal  ****************/ 
receiver  (coded.data,  noise,  limit,  rcvd_tr); 

prlntf  ("Receiver  done  [It ,snr]=EXl,Xl]\n" ,  pe_avg,  snr);  /*  Progress  Indicator  */ 

/t**********  PRINT  OUT  data  for  spectrum  plots  *************/ 

If  (saveData  ==  1) 

■E 

//  For  Sundevll 
If  (where_runnlng  ==  1) 

■E 

lnterferer_out  = 
f  open 

("x:\\documents\\research\\converslon\\results\\lnt erf erer_out.txt" , 
"w"); 

nolse_out  = 
f  open 

("x: \\documents\\research\\converslon\\results\\nolse_out . txt" , 

"w"); 
rcvd_tr_out  = 
f  open 

("x: \\documents\\research\\converslon\\results\\rcvd_tr_out . txt" , 
"w"); 

a_hat_out  = 
f  open 

("x: \\documents\\research\\converslon\\results\\a_hat_out . txt" , 

"w"); 

} 

//  For  L90 

If  (where_runnlng  ==  2) 

■E 


lnterferer_out  = 

f open  ("c:\\temp\\mlke\\results\\lnterferer_out . txt" , 

"w"); 

nolse_out  = 

f open  ("c:\\temp\\mlke\\results\\nolse_out . txt" ,  "w") ; 
rcvd_tr_out  = 

f open  ("c:\\temp\\mlke\\results\\rcvd_tr_out . txt" ,  "w") ; 
a_hat_out  = 

f open  ("c:\\temp\\mlke\\results\\a_hat_out . txt" ,  "w") ; 

} 

for  (1=0;  1  <  Ntot;  1++) 

■E 

f prlntf  (lnterferer_out ,  "Xf \n",  Interferer [1] ) ; 
f prlntf  (nolse_out,  "Xf\n",  noise [1]); 
fprlntf  (rcvd_tr_out,  "Xf\n",  rcvd_tr[l]); 

} 

fclose  (lnterferer_out) ; 
fclose  (nolse_out) ; 
fclose  (rcvd_tr_out) ; 

saveData  =  0;  //  Don’t  save  data  anymore 

} 

/************  Implement  a  linear  receiver  with  truncation 
decision  (chips,  rcvd_tr,  lc_dec); 

prlntf  ("LC:  Decision  done,  [It ,snr]=[Xl,Xl]\n" ,  pe_avg,  snr);  /*  Progress  Indicator  */ 
for  (1  =  0,  num_errs_lc  =  0.0;  1  <  B;  1++)  /*  Find  no.  of  error  for  linear  correlator  */ 

num_errs_lc  =  num_errs_lc  +  (O.S  *  abs  ((lc_dec[l]  -  refEl]))); 

/t************  Implement  LO  Detector  Techniques  **********************/ 

/*  Find  the  appropriate  whltener  */ 

j  =  (Int)  (floor  (Ntot  /  /  Bar))  -  1);  /*  last  array  element  of  Interest  */ 

temp  =  valloc  (NULL,  j);  /*  allocate  for  passing  In  partial  rcvd_tr  */ 

for  (1=0;  1  <  j ;  1++) 

*(temp  +  1)  =  *(rcvd_tr  +1);  /*  Copy  rcvd_tr[]  Into  tempE]  */ 

modcov  (temp,  j,  P,  a_hat,  error);  /*  Find  the  whltener  */ 
vf ree  (temp) ;  /*  clean  up  */ 

/***********  Implement  a  linear  correlator  on  the  whitened  signal  ♦♦♦♦♦♦♦/ 
temp  =  valloc  (NULL,  1);  /*  create  vector  with  one  element  =  1  ♦/ 

tempEO]  =  1; 
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/*  find  whitened  received  signal  */ 

filter  (vscale  (a_hat,  P  +  1,  -1.0),  P  +  1,  temp,  1,  rcvd_tr, 

Ntot,  whitened); 

vf ree  (temp) ;  /*  clean  up  */ 

decision  (chips,  whitened,  wc_dec);  /*  mahe  decisions  */ 

prlntf  ("LC  on  Wh:  Decision  done,  [lt,snr]=[Xl,Xl] \n" ,  pe_avg,  snr) ;  /*  Progress  */ 

for  (1  =  0,  num_errs_wc  =  0.0;  1  <  B;  1++)  /*  Find  no.  of  errors  for  LC  on  wh.  */ 

num_errs_wc  =  num_errs_wc  +  (0.5  *  abs  ((wc_dec[l]  -  ref[l]))); 

/t*************  Implement  ARLO  and  LO  Receivers  ***************/ 

/*  Preparation  */ 

vlnlt  (g_lo,  Ntot,  O.O);  /*  Initialize  to  0.0  */ 

vlnlt  (g_ar,  Ntot,  0.0); 
vlnlt  (g_lw,  Ntot,  0.0); 

Nsamh  =  Bh  *  Nb;  /*  Number  of  samples  In  each  histogram  */ 

temp2  =  valloc  (NULL,  1);  /*  Allocate  for  a_hat=-l  (P=0)  */ 

temp2[0]  =  -1; 

for  (hlstonum  =  0;  hlstonum  <  (dnt)  floor  (B  /  Bh));  hlstonum++) 

-{  /*  Process  Bh*Nb  samples 

at  a  time  for  a  histogram  */ 

vlnlt  (block,  Nsamh,  0.0);  /*  Inlt  block  to  prepare  for  histogram  */ 

for  (1  =  hlstonum  *  Nsamh,  j  =  0; 

1  <  ((hlstonum  +  1)  *  Nsamh  -  1);  1++,  j++) 

♦(block  +  j)  =  *(rcvd_tr  +1);  /*  Apply  the  correct  symbols  */ 
temp  =  valloc  (NULL,  l) ;  /*  Allocate  temp[o]  =  -1  for  LO  case  */ 

temp[o]  =  -1; 

ar_hlsto  (block,  Bh  *  Nb,  temp,  1,  h_lo,  bw_lo,  mln_lo);  /*  LO  (P=0)  a_0=-l  */ 

ar_hlsto  (block,  Bh  *  Nb,  a_hat,  P  +  1,  h_wh,  bw_wh,  mln_wh) ;  /*  AR  methods  */ 
vfree  (temp) ;  /*  clean  up  */ 

for  (symcnt  =  0;  symcnt  <  Bh;  symcnt++)  /*  Apply  nonlinearity  to  each  symbol  */ 

i 


} 


vlnlt  (symbol,  Nb,  0.0);  /*  Look  at  one  symbol  */ 

1  =  hlstonum  *  Nsamh  +  ((symcnt  -  1)  *  Nb) ;  /*Flnd  Indices  of  current  symbol  */ 
for  (j  =  0;  j  <  Nb;  j++) 

■{ 


♦(symbollndx  +  j)  =  1  +  j; 

♦(symbol  +  j)  =  ♦(rcvd_tr  +  1  +  j); 
♦(fllsym  +  j)  =  ♦(whitened  +  1  +  j); 


/♦  Get  current  symbol  ♦/ 

/♦  Get  current  whitened  symbol  ♦/ 


temp  =  valloc  (NULL,  Nb) ;  /♦  To  hold  g_xx[symbollndx]  ♦/ 

temp  =  nip  (symbol,  Nb,  mln_lo,  h_lo,  K  +  1,  bw_lo,  a_hat_negl,  1,  Nb) ;  /♦  lo  detector 


<  Nb;  1++,  j++) 


,  h_wh,  K  ■ 


for  (1  =  symbollndx[0] ,  j  =  0; 

*(g_lo  +  1)  =  ♦(temp  +  j); 
temp  =  nip  (fllsym,  Nb,  mln_ 
for  (1  =  symbollndxCo] »  j  =  0;  j  <  Nb;  1++,  j++) 
*(g_lM  +  1)  =  ♦(temp  +  j); 
temp  =  nip  (symbol,  Nb,  mln_w] 
for  (1  =  symbollndx[0] ,  j  =  0: 

*(g_lM  +  1)  =  ♦(temp  +  j); 


/♦  Assign  g_lo  =  temp  ♦/ 


1,  bw_wh,  a_hat_negl,  1,  Nb) ;  /♦  Iw  detector 
/♦  Assign  g_lw  =  temp  ♦/ 


h_wh,  K  +  1,  bw_wh, 
j  <  Nb;  1++,  j++) 


a_hat,  P  +  1,  Nb) ;  /♦  ar  detector 

/♦  Assign  g_ar  =  temp  ♦/ 


♦/ 


♦/ 


♦/ 


} 

vfree  (temp2) ; 

/♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦  Make  ARLO  and  LO  Detector  Decisions  ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦/ 
decision  (chips,  g_lo,  lo_dec);  /♦  make  decisions  ♦/ 

prlntf  ("LO:  Decision  done,  [It ,snr]=[Xl,Xl]\n" ,  pe_avg,  snr);  /♦  Progress  ♦/ 
decision  (chips,  g_lw,  lw_dec); 

prlntf  ("LW:  Decision  done,  [It ,snr]=[Xl,Xl]\n" ,  pe_avg,  snr);  /♦  Progress  ♦/ 
decision  (chips,  g_ar,  ar_dec); 

prlntf  ("AR:  Decision  done,  [It ,snr]=[Xl,Xl]\n" ,  pe_avg,  snr);  /♦  Progress  ♦/ 

for  (1  =  0,  num_errs_lo  =  0.0,  num_errs_lw  =  0.0,  num_errs_ar  =  0.0;  1  <  B;  1++)  /♦  Count  number  of  errors  ♦/ 

■[ 


num_errs_lo  =  num_errs_lo  +  (0.5  ♦  abs  ((lo_dec[l]  -  ref[l]))) 
num_errs_lw  =  num_errs_lw  +  (0.5  ♦  abs  ((lw_dec[l]  -  ref[l]))) 
num_errs_ar  =  num_errs_ar  +  (0.5  ♦  abs  ((ar_dec[l]  -  ref[l]))) 


/♦  lo  errors  ♦/ 
/♦  Iw  errors  ♦/ 
/♦  ar  errors  ♦/ 


} 

//vf ree(lo_dec) ;vfree(lw_dec) ;vf ree(ar_dec) ; 

/*************  Compute  probability  of  bit  error  for  this  MC  Iteration  and  this  snr  ♦/ 

Pb_lc [pe_avg] [snr]  =  num_errs_lc  /  B;  /♦  linear  correlator  ♦/ 

Pb_wc [pe_avg] [snr]  =  num_errs_wc  /  B; 

Pb_lo[pe_avg] [snr]  =  num_errs_lo  /  B; 

Pb_lw[pe_avg] [snr]  =  num_errs_lw  /  B; 

Pb_ar[pe_avg] [snr]  =  num_errs_ar  /  B; 

prlntf  ("Monte  Carlo  Iteration  Xl  done,  [snr]=[Xl] \n\n",  pe_avg,  snr);  /♦  Progress  Indicator  ♦/ 
}  /♦  End  Monte  Carlo  Iteration  ♦/ 

prlntf  ("SNR  Iteration  Xl  done.  \n\n",  Eb_mln  +  snr);  /♦  Progress  Indicator  ♦/ 

}  /♦  End  SNR  Iteration  ♦/ 

/t*************  Compute  average  P[err]  from  MC  analysis  ♦♦♦♦♦♦♦♦♦♦♦/ 
temp  =  valloc  (NULL,  5);  /♦  5  different  receivers  ♦/ 

for  (1=0;  1  <  (Eb_max  -  Eb_mln  +  1);  1++) 

■[ 


vlnlt  (temp,  5,  0.0);  /♦  Initialize  temp  to  zeros  ♦/ 

for  (j  =  0;  j  <  It;  j++)  /♦  Sum  up  over  all  MC  Iterations  ♦/ 
■[ 


} 


temp[o]  =  temp[o] 
temp[l]  =  temp[l] 
temp [2]  =  temp [2] 
temp [3]  =  temp [3] 
t  emp  [4]  =  t  emp  [4] 


Pb_lc[j] [1] ; 
Pb_wc[j] [1] ; 
Pb_lo[j] [1] ; 
Pb_lw[j] [1] ; 
Pb_ar[j] [1] ; 


/♦  linear  correlator  ♦/ 


uselnput_  =1;  /♦  Perform  In-place  ♦/ 

temp  =  vscale  (temp,  5,  1.0  /  It);  /♦  Average  over  MC  ♦/ 

uselnput_  =0;  /♦  Reset  ♦/ 
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} 


Pb_lc_avg[l]  =  temp[0]; 

Pb_wc_avg[l]  =  temp[l]; 

Pb_lo_avg[l]  =  temp[2]; 

Pb_lw_avg[l]  =  temp[3]; 

Pb_ar_avg[l]  =  temp[4]; 

} 

vfree  (temp) ;  /*  Clean  up  */ 

/t***********  Save  to  disk  P[err]  info  *************/ 

//  For  sundevil 
if  (where_running  ==  1) 

■{ 


snr_vals  = 

f open  ("x : \\documents\\research\\conversion\\results\\snr_vals_ . txt" , 
■'W'); 


Pblc  = 

f open  ("x:\\documents\\research\\conversion\\results\\pblc . txt", 
Pbwc  = 

f open  ("x:\\documents\\research\\conversion\\results\\pbwc . txt", 
Pblo  = 

f open  ("x :  WdocumentsWresearchWconversionWresultsWpblo .  txt" , 
Pblw  = 

f open  ("x :  WdocumentsWresearchWconversionWresultsWpblw .  txt" , 
Pbar  = 

f open  ("x :  WdocumentsWresearchWconversionWresultsWpbar .  txt" , 

} 

//  For  L90 

if  (where_running  ==  2) 

■{ 

snr_vals  =  fopen  ("c :\\temp\ \mike\\results\\snr_ vals_.txt" ,  "w"); 
Pblc  =  fopen  ("c :\\temp\\mike\\results\\pblc .txt",  "w"); 

Pbwc  =  fopen  ("c :\\temp\\mike\\results\\pbwc .txt",  "w"); 

Pblo  =  fopen  ("c iWtempWmikeWresultsWpblo.txt",  "w"); 

Pblw  =  fopen  ("c :\\temp\\mike\\results\\pblw.txt",  "w"); 

Pbar  =  fopen  ("c iWtempWmikeWresultsWpbar.txt",  "w"); 

} 

for  (i  =  0;  i  <  (Eb_max  -  Eb_min  +  1);  i++) 

■{ 


} 


f printf 

(snr_vals,  "Xi\n",  Eb_min  +  1); 

/*  Save  snr  vector  */ 

f printf 

(Pblc , 

"Xf\n" 

Pb_lc_avg[l] ) ; 

/*  linear  correlator  */ 

f printf 

(Pbwc , 

"Xf\n" 

Pb_wc_avg[l] ) ; 

f printf 

(Pblo, 

"Xf\n" 

Pb_lo_avg[l] ) ; 

f printf 

(Pblw, 

"Xf\n" 

Pb_lw_avg[i] ) ; 

f printf 

(Pbar, 

"Xf\n" 

Pb_ar_avg[i] ) ; 

/*  Clean  house  */ 
f close  (Pblc) ; 
f close  (Pbwc) ; 
f close  (Pblo) ; 
f close  (Pblw) ; 
f close  (Pbar) ; 
f close  (snr_vals) ; 
mathfree  0 ; 


"w") 

"w") 

"w") 

"w") 

"w") 


A. 2. 2  AR  Histo.c 


/*  AR_histo.c 


by  Mike  Banys 
March  3,  2001 


*  This  estimates  the  pdf  of  a  white  noise  sequence  that  drives  a  P-th 

*  order  auto-regressive  filter  model.  It  returns  the  vector  function 

*  hO  which  is  not  the  pdf,  but  rather  the  natural  log  of  the 

*  derivative  of  the  pdf,  to  be  used  to  the  LO  non-linear  processor. 


*/ 

#include  <stdlib.h> 

#include  <stdio.h> 

#include  <math.h> 

#include  "mathlib . h" 

#include  "slm_parameters .h" 

#include  "sim_parameters_extern.h" 

voidar_histo  (Real_Vector  stimulus,  unsigned  stimulus_length,  Real_Vector  a, 
unsigned  a.length, 

Real_Vector  h,  Real_Vector  bin_width,  Real_Vector  mini) 

/*  bin_wldth  ft  mini  are  vectors  to  enable  pass-by-reference  */ 

■{ 

int  i,  j;  /*  counter  variable  */ 

unsigned  N;  /*  samples  for  histogram  */ 

Real_Vector  white,  f ilter_a_variable,  /*  used  for  sending  -1  into  filterO  */ 

pdf,  /*  pdf  of  white  noise  (frequency  count  of  histogram)  */ 

bp,  /*  histogram  bin  centers  */ 

neg_a;  /*  Temporary  variable  to  hold  -a  */ 

white  =  valloc  (NULL,  stimulus_length) ; 
if  (a_length  ==  1) 

vcopy  (white,  stimulus,  stimulus_length) ;  /*  using  ARLO  */ 
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/*  using  filtered  LO  */ 

■{ 

f ilter_a_variable  =  valloc  (NULL,  1); 

f ilter_a_variable[0]  =  1; 

neg_a  =  vscale  (a,  a_length,  -1.0); 

filter  (neg_a,  a_length,  f liter. a_variable,  1,  stimulus, 
stimulus.length,  white) ; 

vfree  (neg_a) ; 

} 

/*  Build  the  pdf  estimate  via  the  histogram  with  K  bins  */ 

N  =  Nb  *  Bh;  /*  samples  calculated  for  this  histogram  */ 

pdf  =  valloc  (NULL,  K) ; 
bp  =  valloc  (NULL,  K) ; 

hist  (white,  stlmulus.length,  K,  pdf,  bp);  /*  Compute  histogram  estimate  of  white  noise  */ 
for  (i  =  0;  1  <  K;  i++)  /♦  Find  possible  negative  values  or  zeros  */ 

■{ 

if  (*(pdf  +  i)  <=  0.0)  /*  set  them  equal  to  1.0  for  calculation  of  g  */ 

♦(pdf  +  i)  =  1.0; 

} 

bin_width[o]  =  (Real)  ♦  (bp  +  1)  -  ♦bp;  /♦  Width  of  each  equi-spaced  bin,  typecasted  to  Real  ♦/ 

mini[o]  =  vminval  (white,  stimulus.length,  &i);  /♦  Minimum  value  in  the  histogram  ♦/ 

useinput.  =1;  /♦  overwrite  pdf[]  with  scaled  pdf[]  ♦/ 

pdf  =  vscale  (pdf,  K,  (1  /  (♦bin.width  ♦  N)));  /♦  scale  histogram  so  that  pdf  integrates  to  1  ♦/ 

useinput.  =0;  /♦  reset  ♦/ 

/♦  The  "h"  function  Is  the  derivative  of  the  In  of  the  pdf  of  the  white  noise  sequence  ♦/ 

♦h  =  (l.O  /  ♦bin.width)  ♦  (log  (♦(pdf  +  1))  -  log  (♦pdf)); 
for  (i=l;i<K-l;  i++) 

♦(h  +  i)  = 

(O.S  /  ♦bin.width)  ♦  (log  (♦(pdf  +  1  +  1))  -  log  (♦(pdf  +  1  -  1))); 

♦(h  +  K  -  1)  = 

(l.O  /  ♦bin.width)  ♦  (log  (♦(pdf  +  K  -  1))  -  log  (♦(pdf  +  K  -  2))); 

♦(h  +  K)  =  0.0;  /♦  For  received  values  outside  support  of  h  (a  junkbin)  ♦/ 

/♦  Clean  House  ♦/ 

vfree  (f ilter_a_variable) ; 
vfree  (pdf) ; 
vfree  (bp) ; 
vfree  (white) ; 


A. 2. 3  AR_Theory.c 


/♦AR  Theory 

♦  by  Fernando  Martinez  Vallina 

* 

♦  Calculates  the  pdf  of  a  white  noise  sequence  that  drives  the  F-th 

♦  auto-regresslve  filter  model.  Returns  a  function  hO  that  Is  the 

♦  natural  log  of  the  derivative  of  the  pdf. 


♦  Inputs;  K  -  Number  of  data  points 

♦  T  -  Support  of  the  pdf 

♦  iidtype  -  Type  of  the  white  noise  pdf 

♦  lldparms  -  Parameters  of  the  pdf 

♦  ISR  -  Interferer  to  Signal  Ratio  (dB) 

♦ 

♦  Outputs:  h  -  d/d(rho)  of  ln(fw()) 


♦  This  Is  based  on  the  original  AR  Theory  code  In  MATLAB 
♦/ 

#include  <math.h> 

#include  <stdlib.h> 

#include  "mathlib . h" 

#include  "slm.parameters .h" 

#include  "sim_parameters_extern.h" 
void 

ar.theory  (int  iidparmsH,  char  iidtype[],  int  ISR,  Real  K, 

Real  T,  Real.Vector  h) 

■{ 

Real  xstep,  xmax,  slgl,  power; 

Int  length,  1;  //length  Is  the  length  of  the  sequence,  1  Is  a  counting  variable 

Real.Vector  x,  mu,  sig; 

power  =  pow  (10,  (ISR  /  10)); 

xmax  =  T  /  2; 

xstep  =  T  /  (K  -  1) ; 

length  =  T  /  xstep; 

X  =  valloc  (NULL,  length); 
x[o]  =  -xmax; 

for  (i  =  1;  i  <=  length;  i++) 

■{ 

x[i]  =  x[i  -  l]  +  xstep; 

} 

h  =  valloc  (NULL,  length); 
if  (iidtype  ==  ’gs’) 

■{ 


mu  =  valloc  (NULL,  length) ; 
for  (1=0;  1  <  length;  i++) 
■{ 
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muCi]  =  iidparms[l]; 


} 

slg  =  valloc  (NULL,  length); 
for  (1=0;  1  <  length;  i++) 

■{ 

slg[l]  =  lldpams  [2] ; 

} 

for  (1=0;  1  <  length;  1++) 

■{ 

h[l]  =  -(x[l]  -  muCl])  /  (power  ♦  slg[l]  ♦  slg[l]); 

} 

} 

If  (lldtype  ==  ’Ip’) 

■{ 

slgl  =  lldparns [1] ; 

slgl  =  slgl  *  slgl  *  power; 

} 

vfree  (x) ; 
vfree  (mu) ; 
vfree  (slg) ; 

} 


A. 2. 4  Contwave.c 

/*  Interferer  must  be  Initialized  to  zero  (or  contain  values  from 
another  Interference  type)  before  calling  this  function 
*/ 

#lnclude  <math.h> 

#lnclude  <stdllb.h> 

#lnclude  <stdlo.h> 

#lnclude  "mathllb . h" 

#lnclude  "slm_parameters_extern.h" 
void 

contwave  (Real_Vector  Interferer) 

■{ 

Real  amp,  /*  Individual  Interferer  amplitude  */ 

delta_t;  /*  time  Interval  */ 

Int  1,  Intf;  /*  counter  */ 

Real_Vector  time,  /*  Time  sequence  */ 

phase;  /*  phase  of  each  Interferer  */ 

double  pi  =  3.1415926535; 

amp  =  sqrt  (2  ♦  pow  (10,  ISR  /  10)  /I);  /*  Individual  Interferer  amplitudes  */ 

phase  =  valloc  (NULL,  I); 

for  (1=0;  1  <  I;  1++)  /♦  Set  each  Interferer’s  phase  to  0  */ 

♦(phase  +1)  =0; 

delta_t  =  (Real)  1.0  /  Ntot;  /*  Time  Interval  */ 
time  =  valloc  (NULL,  Ntot);  /*  Time  sequence  */ 
for  (1=0;  1  <  Ntot;  1++) 

♦(time  +  1)  =  delta_t  ♦  1; 

for  (Intf  =  0;  Intf  <  I;  lntf++)  /♦  Create  sinusoid  Interferer  ♦/ 

■{ 

for  (1=0;  1  <  Ntot;  1++) 

♦(Interferer  +  1)  = 

♦(Interferer  +  1)  + 

(amp  ♦  sin  (2  ♦  pi  ♦  freqClntf]  ♦  timed]  +  phaseClntf] )) ; 

} 

vfree  (time) ; 
vfree  (phase) ; 


A. 2. 5  Conv.c 


/♦  Convolution  Function 
* 

♦  by  Michael  Banys 

* 

♦  conv(vl,  lengthl,  v2,  length2,  output) 

* 

♦  Inputs:  vl  =  Vector  1 

♦  v2  =  Vector  2 

♦  lengthl  =  length(vl) 

♦  length2  =  length(v2) 

♦ 

♦  Outputs:  output  =  convolution  of  vecl  with  vec2 

* 

♦  "The  resulting  vector  Is  output_length=lengthl+length2-l . 

♦  If  vl  and  v2  are  vectors  of  polynomial  coefficients,  convolving 

♦  them  Is  equivalent  to  multiplying  the  two  polynomials." 

♦ 

♦  This  Is  based  on  MATLAB’s  conv  function. 

♦/ 

#lnclude  <stdllb.h> 

#lnclude  "mathllb . h" 
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#include  "sini_paraiaeters_extern.h" 

voidconv  (Real_Vector  vl,  Int  lengthl,  Real_Vector  v2, 

int  length2,  Real_Vector  output) 


/*  Use  same  format  as  fllterO  for  this  operation. 

*  conv(a,b)  Is  the  same  as  conv(b,a),  but  we  can  make  It 

*  go  substantially  faster  If  we  swap  arguments  to  make  the  first 

*  argument  the  shorter  of  the  two. 

*/ 

Int  1,  /*  Counter  variable  */ 

f llter_order,  /*  Order  of  numerator  */ 

f llter_lteratlon,  /*  Variable  tracking  loop  through  Input  data  */ 


element.num, 

f lrst_length, 
second.length; 
Real_Vector  v_alt, 
*f lrst_v, 

*second_v, 
f llter_delay; 
Real  sumx, 


If  (lengthl  >  length2) 
■{ 


/*  Variable  Indicating  current  coefficient  and  delay 
being  processed  */ 

/*  Length  of  vector  that  Is  first  Into  conv  */ 

/*  Length  of  vector  that  Is  used  as  data_ln  */ 

/*  Temp  vector  for  zero-padding  vl  or  v2  */ 

/*  Pointer  for  first  argument  Into  convolution, 

as  explained  below.  Nay  be  assigned  to  vl  oi 

/*  as  above  */ 

/*  Filter  delay  values  */ 

/*  Temporary  variable  representing  sum  of  elements  of  a[] 
multiplied  by  elements  of  filter. delay[]  */ 

/*  Temporary  variable  reperesentlng  sum  of  elements  of  b[] 
multiplied  by  elements  of  f llter_delay[]  */ 


If  (length2  >  l) 

■{ 

v_alt  =  valloc  (NULL,  (lengthl  +  length2  -  1));  /*  For  zero-padding  vl  */ 

for  (1=0;  1  <  lengthl;  1++) 

*(v_alt  +  1)  =  *(vl  +1);  /*  Copy  vl  Into  v_alt  */ 

for  (1;  1  <  (lengthl  +  length2  -  1);  1++)  /♦  zero  pad  rest  of  v_alt  */ 

*(v_alt  +1)  =0; 


flrst.v  =  &v2;  /*  v2  will  be  first  vector  Into  the  conv  operation  */ 

second.v  =  &v_alt; 
flrst.length  =  length2; 
second.length  =  lengthl; 


} 


■{ 

If  (lengthl  >  l) 

■{ 

v_alt  =  valloc  (NULL,  (lengthl  +  length2  -  1));  /*  For  zero-padding  v2  */ 

for  (1=0;  1  <  length2;  1++) 

*(v_alt  +  1)  =  ♦(v2  +1);  /*  Copy  v2  Into  v_alt  */ 

for  (1;  1  <  (lengthl  +  length2  -  1);  1++) 

*(v_alt  +  1)  =  0;  /*  Zero  pad  rest  of  v_alt  */ 

} 

flrst.v  =  &vl;  /*  vl  will  be  first  vector  Into  the  conv  operation  */ 

second.v  =  &v_alt; 
flrst.length  =  lengthl; 
second.length  =  length2; 

} 

/*  Begin  actual  convolution  using  Direct  Form  II  structure,  as  In  fllterO  */ 
fllter.delay  =  valloc  (NULL,  f Irst.length) ;  /*  Allocate  memory  for  fllter.delay  length  */ 

fllter.order  =  flrst.length;  /*  Set  order  of  numerator  */ 

for  (1=0;  1  <  f Irst.length;  1++)  /*  Set  w[n]  delay  array,  filter .delay,  to  all  zeros  */ 

*(fllter.delay  +1)  =0; 

for  (f llter.lteratlon  =  0;  f llter.lteratlon  <  (first. length  +  second. length  -  1);  f liter. lteratlon++)  /*  Start  Filtering  */ 

i 

sumx  =  sumy  =  0.0;  /*  Initialize  variables  */ 

/*  Can  SKIP  summing  left  side,  as  there  are  no  denominator  coefficients  */ 

/*  Sum  right  side  of  Direct  Form  II  Structure  */ 

for  (element. num  =  1;  element. num  <  first. length;  element.num++) 

sumy  =  sumy  +  (*(*flrst.v  +  element. num)  *  f liter. delay [element .num] ) ; 

/*  Compute  current  filter  delay  (filter. delayEo] )  */ 

♦filter. delay  =  *(*second.v  +  f llter.lteratlon)  +  sumx; 

/*  Compute  output  element,  y(n)  (y[f llter.lteratlon] )  */ 
output [filter. Iteration]  =  *(*flrst.v)  ♦  filter .delay[0]  +  sumy; 

/*  Delay  each  filter. delay[]  element  by  one  */ 

for  (element. num  =  first. length  -  1;  element. num  >  0;  element. num — ) 
f liter .delayCelement. num]  =  filter. delayCelement.num  -  1] ; 

} 

vfree  (f llter.delay) ; 
vfree  (v.alt) ; 


A. 2. 6  Decision. c 


#lnclude 

#lnclude 

#lnclude 

#lnclude 

void 

decision 


<stdllb.h> 

<math.h> 

"mathllb.h" 

"slm.parameters.extern.h" 

(Real. Vector  pn.  Real. Vector  g.  Real. Vector  dec) 


v2  ♦/ 
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int  i,  j ,  k; 

Real  corr; 

for  (1  =  0,  k  =  0;  1  <  Ntot;  1=1+  Nt,  k++) 


/*  Counter  */ 

/*  Correlation  of  pn  sequence  with  g  */ 


■{ 

for  (j  =  0,  corr  =  0.0;  j  <  Nh;  j++)  /*  Correlate  each  symbol  with  pn  and  SUM  each  correlation  */ 

corr  =  corr  +  (pn[l  +  j]  *  g[l  +  j]); 

If  (corr  <=  0) 

♦(dec  +  k)  =  -1; 


} 


} 


else 
♦  (dec 


k)  =  l; 


A.2.7  DSSS.c 


#lnclude  <tlme.h> 

#lnclude  <stdllb.h> 

#lnclude  "mathllb . h" 

#lnclude  "slm_parameters_extern.h" 

volddsss  (Real_Vector  ref,  Int  amp,  Real_Vector  coded. data,  Real.Vector  chips) 

■{ 

unsigned  Int  1,  /♦  temp  var  ♦/ 

j ,  /♦  temp  var  ♦/ 

Real.Vector  data,  tmp; 

long  Itlme;  /♦  For  random  seed  ♦/ 

data  =  valloc  (NULL,  Ntot);  /♦  Sample  each  symbol  Nb  times  (Nb  samples  per  symbol)  ♦/ 

for  (1  =  0,  k  =  0;  1  <  Ntot;  1=1+  Nb,  k++)  /♦  compute  data  vector  ♦/ 

■{ 

for  (j  =  0;  j  <  Nb;  j++) 

♦(data  +  1  +  j)  =  ♦(ref  +  k) ; 

} 

tmp  =  valloc  (NULL,  B  ♦  PG) ; 

Itlme  =  time  (NULL) ; 

srand  ((unsigned)  Itlme  /  2); 

for  (1=0;  1  <  B  ♦  PG;  1++)  /♦  Create  PN  chipping  sequence  ♦/ 

■{ 

♦(tmp  +  1)  =  (2  ♦  (Real)  rand  0  /  RAND.MAX)  -  1.0;  /♦  PN  sequence  B^PG  chips  long  ♦/ 

If  (♦(tmp  +  1)  <  0)  /♦  Random  chips  at  1  sample  per  chip  ♦/ 

♦  (tmp  +  1)  =  -1.0; 
else  If  (♦(tmp  +1)  >0) 

♦  (tmp  +  1)  =  1.0; 
else 

♦ (tmp  +1)  =0.0; 

} 

for  (1=0, k=0;  1<  Ntot;  1  =  1  +  Nc,  k++)  /♦  Sample  each  chip  Nc  times  ♦/ 

■{ 

for  (j  =  0;  j  <  Nc;  j++) 

♦(chips  +  1  +  j)  =  ♦(tmp  +  k) ; 

} 

for  (1=0;  1  <  Ntot;  1++)  /♦  The  sampled  signal  modulated  by  the  PN  sequence  ♦/ 

♦(coded. data  +  1)  =  amp  ♦  (♦(chips  +  1))  ♦  (♦(data  +1)); 
vfree  (data) ; 
vfree  (tmp) ; 


A.2.8  FFTShift.c 


/♦  FFTSHIFT.C 
* 

♦  Swaps  the  left  and  right  halves  of  the  Input  vector,  placing  the  DC 

♦  component  In  the  middle  of  the  spectrum. 

♦  l.e.  a  =  [0  1  2  3  4]  — >  [3  4  0  1  2] 

♦  l.e.  a  =  [0  1  2  3  4  5]  — >  [  3  4  5  0  1  2] 

♦ 

♦  vl  =  Input  complex  vector 

♦  lengthl  =  length  of  complex  vector 


♦  by  Michael  Banys 

♦  February  7,  2000 
♦/ 

#lnclude  <stdllb.h> 

#lnclude  <math.h> 

#lnclude  "mathllb . h" 

#lnclude  "slm.parameters.extern.h" 

voldfftshlft  (Complex.Vector  vl,  unsigned  Int  lengthl) 

■[ 

unsigned  Int  k, 
j,  midpoint; 

Complex  tempi, 
temp2 ; 

midpoint  =  (unsigned  Int)  cell  ((double)  lengthl  /  2) 
/♦  If  lengthl  Is  odd,  this  finds  exact  midpoint. 


/♦  counter  ♦/ 

/♦  Middle  of  vl  ♦/ 

/♦  temp  Storage  ♦/ 

-  1;  /♦  Find  midpoint  of  vl  ♦/ 
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*  If  lengthl  Is  even,  this  find  the  lesser  of  the  two  midpoints  */ 

If  (fmod  ((double)  lengthl,  2.0))  /♦If  vl  contains  ODD  number  of  elements  ♦/ 

■{ 


} 


temp2.r  =  vl [midpoint] .r; 
temp2.1  =  vl [midpoint] . 1; 

for  (k  =  midpoint,  j  =  (lengthl  -  1);  k  >=  1;  k — ,  j — )  /♦  Iterate  from  midpoint  to  zero  ♦/ 
■{ 


tempi .r  =  vl  [j] .r; 
tempi . 1  =  vl  [j] .1; 
vl  [j] .r  =  temp2.r; 
vl  [j] . 1  =  temp2. 1; 
temp2.r  =  vl[k  -  l].r; 
temp2.1  =  vl[k  -  l].l; 
vl[k  -  l] .r  =  tempi. r; 
vl[k  -  1] .1  =  tempi. 1; 

} 

vl [midpoint] .r  =  temp2.r; 
vl [midpoint] . 1  =  temp2.1; 


/♦  Save  first  storage  location  ♦/ 
/♦  Move  a  value  ♦/ 

/♦  Save  second  storage  location  ♦/ 
/♦  Move  a  value  ♦/ 


} 


/♦  If  vl  contains  EVEN  number  of  elements  ♦/ 
■{ 


for  (k  =  midpoint,  j  =  (lengthl  -  1);  k  >  0;  k — ,  j — ) 


} 


tempi .r  =  vl  [j] .r 
tempi .1  =  vl  [j] . 1 
vl [j] .r  =  vl  [k] .r 
vl [j] . 1  =  vl  [k] . 1 
vl [k] .r  =  tempi .r 
vl [k] . 1  =  tempi . 1 


tempi .r  =  vl [j] 
tempi . 1  =  vl [j] 
vl[j]  .r  =  vl[k] 
vl[j]  .1  =  vl[k] 
vl  [k] .r  =  tempi .r; 
vl  [k] . 1  =  tempi . 1 


/♦  Save  value  ♦/ 
/♦  Store  value  ♦/ 
/♦  Store  value  ♦/ 

/♦  Save  value  ♦/ 
/♦  Store  value  ♦/ 
/♦  Store  value  ♦/ 


A. 2. 9  Filter. c 


/♦  filter. c 

♦  Version  2.0,  5/1/01 

♦  Written  by  Michael  Banys  &  Angela  Kaczmarskl 

♦ 

♦  Revision  Control: 

♦  2.0  Changed  to  allow  a  denominator  of  1 

♦  1.0  Original  version 

♦ 

♦/ 

/♦  OUTPUT  =  fllter(  b[M],  M,  a[N] ,  N,  data_ln[lnput_length] ,  lnput_length,  y[]  ) 

♦ 

♦  b[M]  =  coefficients  of  numerator  of  Impulse  response 

♦  a[N]  =  coefficients  of  denominator  of  Impulse  response 

♦  =  order  of  filter 

♦ 

♦  This  function  computes  the  convolution 

♦  of  vector  "data_ln"  with  the  causal  filter 

♦  whose  time-domain  filter  coefficients 

♦  are  of  the  following  format: 

♦ 

♦  a(l)^y(n)  =  -SUM(k=l,N)  [  a(k)^y(n-k)  ] 

♦  +SUM(k=0,M)  [  b(k)^x(n-k)  ] 

♦ 

♦  If  a(0)  !=  1,  FILTER  will  normalize  the  filter 

♦  coefficients  by  a(0) .  If  a(0)=0,  return  ERROR. 

♦ 

*/ 

#lnclude  <stdllb.h> 

#lnclude  "mathllb . h" 

#deflne  ORDER.ERROR  "Filter  numerator  order  Is  greater  than  denominator  orderXn" 
voldfllter  (Real_Vector  b,  Int  M, 


{ 

Int  f llter_order, 

f llter_lteratlon  =  0, 
element.num, 

1  =  0; 

Real  normallzatlon.value, 
sumx, 

sumy; 

Real_Vector  f llter_delay; 
If  (N  >=  M) 


Real_Vector  a,  Int  N, 

Real_Vector  data_ln,  Int  lnput_length,  Real_Vector  y) 


/♦  Order  of  acting  filter  ♦/ 

/♦  Variable  tracking  loop  through  Input  data  ♦/ 

/♦  Variable  Indicating  current  coefficient  and  delay  being 
processed  ♦/ 

/♦  Iteration  variable  ♦/ 


/♦  a[0]  coefficient  value  ♦/ 

/♦  Temporary 
/♦  Temporary 


/♦  Filter  delay  values  ♦/ 

//  added 


variable  representing  sum  of  elements  of  a[] 

multiplied  by  elements  of  filter _delay[] 
variable  reperesentlng  sum  of  elements  of  b[] 
multiplied  by  elements  of  filter _delay[] 


♦/ 

♦/ 


58 


■{ 

filter_delay  =  valloc  (NULL,  N) ;  /*  Allocate  memory  for  filter_delay  length  */ 

filter_order  =  N;  /*  Set  order  of  filter  */ 

} 

else  //  added 

■{ 

filter_delay  =  valloc  (NULL,  M) ;  /*  Allocate  memory  for  filter_delay  length  *///added 

filter_order  =  M;  /*  Set  order  of  filter  *///  added 

} 

/*  Check  to  see  if  the  length  of  b[]  <  length  of  a[] 

If  length  b[]  >=  length  a[],  output  error  message  and  exit  function  */ 
if  ((M  >=  N)  &  (N  !=  1))  //  added  (N!=l) 

printf  (ORDER.ERROR) ; 

/*  Check  to  see  if  denominator  of  IR  is  normalized  (i.e.  does  a[o]=l 
If  denominator  is  not  normalized,  normalize  IR  by  a[0]  */ 
if  (a[0]  !=  1) 

■{ 

normalization_value  =  a[0] ; 

/*  Normalize  denominator  */ 
for  (i  =  0;  i  <  N;  i++) 

a[i]  =  a[i]  /  normalization_value; 

/*  Normalize  numerator  */ 
for  (i  =  0;  i  <  M;  i++) 

b[i]  =  b[i]  /  normalization_value; 

} 

/*  Set  w[n]  delay  array,  filter_delay[] ,  to  all  zeroes  */ 
for  (i  =  0;  i  <  f ilter_order;  i++)  //  changed  N  to  filter_order 

f ilter_delay[i]  =  0; 

/*  Start  filtering  process  */ 

for  (f ilter_iteration  =  0;  f ilter_iteration  <  input_length; 
filter. it eration++) 

■{ 

sumx  =  sumy  =  0.0;  /*  Initialize  variables  */ 

/*  Sum  left  side  of  Direct  Form  II  Structure  */ 
for  (element.num  =  1;  element.num  <  N;  element_num++) 

sumx  =  sumx  -  a[element_num]  *  f ilter_delay[element_num] ; 

/*  Sum  right  side  of  Direct  Form  II  Structure  */ 
for  (element.num  =  1;  element.num  <  M;  element_num++) 

sumy  =  sumy  +  b[element_num]  *  f ilter_delay[element_num] ; 

/*  Compute  current  filter  delay  (f ilter_delay[0] )  */ 
filter .delay [0]  =  data.inCfilter.iteration]  +  sumx; 

/*  Compute  output  element,  y(n)  (y[f ilter.iteration] )  */ 
yCfilter. iteration]  =  b[0]  *  f ilter.delayCo]  +  sumy; 

/*  Delay  each  f ilter.delay[]  element  by  one  */ 

for  (element. num  =  filter. order  -  1;  element. num  >  0;  element. num — ) 
f ilter.delayCelement.num]  =  f ilter.delayCelement.num  -  l] ; 

} 

vfree  (f ilter.delay) ; 


A. 2. 10  FlipLR.c 

/♦♦♦This  is  the  C  version  of  the  Matlab  FliplrO  function^^^^ 

♦♦♦written  by  Tayo  Ihimoyan,  04/2000^^^^^^/ 

#include  "cmt .h" 

#include  "mathlib . h" 

/♦ - ♦/ 

Real.Matrix 

fliplr  (Real.Matrix  M,  unsigned  m,  unsigned  n) 

■{ 

/♦ 


fliplrO  returns  a  m-by-n  matrix  M  whose  elements 

are  taken  from  M(m-by-n)  with  columns  flipped  in  the  left-right  direction, 
that  is,  about  a  vertical  axis  . 

On  exit,  the  output  matrix  M[] []  is  the  value  of  this  routine. 

♦/ 

/♦ - ♦/ 

unsigned  1,  j,  k; 

Real  1; 

if  (M  ==  NULL) 

■{ 

matherr.  ("fliplr",  E.NULLPTR) ; 
return  NULL; 

} 

if  (n  X  2  ==  0) 
k  =  n  /  2; 

else 

k  =  (n  -  1)  /  2; 

1  =  0.0; 

for  (i  =  0;  i  <  k;  1++) 

■{ 

for  (j  =  0;  j  <  m;  j++) 

■{ 

1  =  M[j][i]; 

M[j][i]  =  M[j][n  -  1  -  1]; 

M[j]  [n  -  1  -  1]  =  1; 
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} 


} 

return  M; 


} 


A. 2. 11  FlipUD.c 

/♦♦♦This  Is  the  C  version  of  the  Matlab  FlipudO  function^^^^ 

♦♦♦written  by  Tayo  Ihimoyan,  04/2000^^^^^^/ 

#include  "cmt .h" 

#include  "mathlib . h" 

/♦ - ♦/ 

Real.Natrix  flipud  (Real.Matrix  M,  unsigned  n,  unsigned  n) 

■{ 

/♦ 


flipudO  returns  a  n-by-n  matrix  M  whose  elements 

are  taken  from  M(m-by-n)  with  rows  flipped  in  the  up-down  direction, 
that  is,  about  a  horizontal  axis  . 

On  exit,  the  output  matrix  M[] []  is  the  value  of  this  routine. 

♦/ 

/♦ - ♦/ 

unsigned  1,  j,  k; 

Real  1; 

if  (M  ==  NULL) 

■{ 

matherr.  ("flipud",  E.NULLPTR) ; 
return  NULL; 

} 

if  (m  X  2  ==  0) 

k  =  m  /  2; 

else 

k  =  (m  -  1)  /  2; 

1  =  0; 

for  (1  =  0;  i  <  k;  i++) 

■{ 

for  (j  =  0;  j  <  n;  j++) 

■{ 

1  =  M[i][j]; 

M[i]  [j]  =  M[m  -  1  -  i]  [j] ; 

M[m  -  1  -  1]  [j]  =  1; 

} 

} 

return  M; 

} 


A. 2. 12  Hist.c 


/♦  hist_mike.c 

♦  Version  1.0,  6/29/2001 

♦  Written  by  Michael  Banys 

♦  Revision  Control: 


♦  hist(  Real_Vector  signal,  int  nbins,  Real_Vector  freq,  Real_Vector  bincenter  ) 


f  req 

bincenter 


=  vector  containing  frequency  distribution  of  each  bin 
=  vector  containing  center  location  of  each  bin 


signal 

nbins 


=  signal  (vector)  from  which  to  compute  histogram 
=  number  of  bins  necessary  for  histogram 


♦  The  histogram  function  bins  the  elements  of  signal  Into  the  number  of  bins  specified  by 

♦  nbins.  The  function  returns  the  number  of  elements  In  each  bln  and  the  center  point 

♦  of  each  bln. 

♦/ 

#include  "mathlib . h" 

#include  <stdio.h> 

#include  <stdlib.h> 

voldhlst  (Real.Vector  signal,  unsigned  slgnal.length, 

unsigned  nbins,  Real_Vector  freq,  Real_Vector  bincenter) 


int  j , 
binfound, 
currentbin; 

Real  mln.slgnal,  max.slgnal,  blnwldth; 
unsigned  temp; 

min_signal  =  vminval  (signal,  signal_length,  fttemp) ;  /♦  Find  min(signal)  ♦/ 


/♦  Counter  variables  ♦/ 

/♦  0  if  data  not  distributed  into  a  bin;  1  if 
/♦  current  bin  being  examined  ♦/ 


it  Is  ♦/ 


60 


max.slgnal  =  vnaxval  (signal,  slgnal.length,  fttenp) ;  /*  Find  niax(slgnal)  */ 

blnwldth  =  (max_slgnal  -  mln.slgnal)  /  nblns;  /*  Calculate  bln  width  */ 
blncenter[0]  =  mln.slgnal  +  0.499999  *  blnwldth;  /*  Center  of  first  bln  */ 

/*  use  0.499999  to  properly  sort  mlnCslgnal)  Into  first  bln  */ 
for  (j  =  1;  j  <  nblns;  j++) 

■{ 

*(blncenter  +  j)  =  *(blncenter  +  j  -  1)  +  blnwldth;  /*  Find  centers  of  rest  of  bins  */ 

*(freq  +  j)  =  0.0;  /♦  Initialize  frequency  vector  to  zero  */ 

} 

freq[o]  =  0.0;  /*  complete  Initialization  of  frequency  vector  to  zero  */ 

for  (j  =  0;  j  <  slgnal_length;  j++)  /♦  Distribute  slgnal[]  Into  bins  */ 

■{ 

blnfound  =  0;  /*  reset  blnfound?  variable  */ 

currentbln  =  0; 

If  (*(slgnal  +  j)  >  (blncenterCnblns  -  1]  +  0.5  *  blnwldth))  /*  If  max(slgnal)  falls  outside  of  bln  */ 
■{ 

freqCnblns  -  1]  =  freqCnblns  -  1]  +  1; 
blnfound  =  1; 

} 

while  (blnfound  =  0) 

■{ 

If  ((*(slgnal  +  j)  >  (blncenterCcurrentbln]  -  0.5  *  blnwldth))  && 

(♦(signal  +  j)  <=  (blncenterCcurrentbln]  +  0.5  *  blnwldth))) 

■{ 

freqCcurrentbln]  =  freqCcurrentbln]  +1;  /*  place  signal  data  Into  current  bln  */ 

blnfound  =  1; 

} 

currentbln++; 

} 

} 

freqCnblns  -  l]  =  freqCnblns  -  l]  +  1;  /*  Adjusted  first  blncenter  causes  max(slgnal) 

not  to  be  sorted  Into  last  bln.  So  add  this 
signal  value  Into  last  bln  manually.  */ 


A. 2. 13  Modcov.c 

/*  Inputs: 

*  N  =  length (x) 

*  p  =  model  order 

*  ARparm  =  AR  parameter  vector  to  return  to  caller 

*  error  =  return  to  caller 

♦/ 

#lnclude  "mathllb . h" 

#lnclude  <stdllb.h> 

#lnclude  "slm_parameters .h" 

#lnclude  "slm_parameters_extern.h" 

voldmodcov  (Real.Vector  x,  unsigned  N,  unsigned  p, 

Real_Vector  ARparm,  Real  error) 

■{ 

Int  1,  j ; 

Real_Matrlx  X,  /*  Toeplltz  matrix  of  x  */ 

R,  /*  temp  storage  matrix  for  fllpud/lr  */ 

Rl,  R2; 

Real_Vector  bl,  /*  temp  storage  vectors  */ 

b2,  b_temp; 

bl  =  valloc  (NULL,  (N  -  p)); 
b2  =  valloc  (NULL,  p  +  1) ; 

for  (1  =  p,  j  =  p;  1  <  N;  1++,  j — )  /*  Copy  vectors  for  use  In  Toeplltz  */ 

< 

♦(bl  +  (1  -  p))  =  x[l] ; 

If  (j  >=  0) 

♦  (b2  +  (1  -  p))  =  x[j] ; 

} 

X  =  toeplltz  (bl,  b2,  (N  -  p) ,  (p  +  1)); 

R  =  mxmull  (X,  X,  p  +  1,  N  -  p,  p  +  1) ;  /♦  Generate  the  matrix  of  normal  equations  ♦/ 

vfree  (bl) ; 

vfree  (b2) ;  /♦  Clean  house  ♦/ 

Rl  =  mxalloc  (NULL,  p,  p) ;  /♦  forward  predictor-error  matrix  ♦/ 

for  (l=l;l<p+l;  1++) 

for  (j=i;j<p+i;  j++) 

Rl[l  -  l][j  -  1]  =  R[l][j]; 

R2  =  mxalloc  (NULL,  p,  p) ;  /♦  Backward  predictor-error  matrix  ♦/ 

for  (1  =  0;  1  <  p;  1++) 

for  (j  =  0;  j  <  p;  j++) 

R2[l]  [j]  =  R[l]  [j]  ; 

R2  =  fllplr  (R2,  p,  p) ; 

R2  =  fllpud  (R2,  p,  p) ; 

bl  =  valloc  (NULL,  p) ;  /♦  Forware  error  coefficients  ♦/ 

for  (l=l;l<p+l;  1++) 

♦  (bl  +  1  -  1)  =  R[l]  [O] ; 

b2  =  valloc  (NULL,  p) ;  /♦  Backward  predictor-error  coefficients  ♦/ 

for  (1  =  0;  1  <  p;  1++) 

■{ 

♦  (b2  +  1)  =  R[p  -  (1  +  l)][p]; 

} 

ARparm[0]  =  -1;  /♦  AR  Coefficients  ♦/ 

b_temp  =  vmxmul  (mxlnv  (mxadd  (Rl,  R2,  p,  p) ,  p) ,  vadd  (bl,  b2,  p) ,  p,  p) ; 
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} 


for  (i=l;i<p+l;  i++) 
ARparmCi]  =  *(b_temp 
inxfree  (Rl) ; 
inxfree  (R2) ; 
vfree  (b_temp) ; 
vfree  (bl) ; 
vfree  (b2) ; 

bl  =  valloc  (NULL,  p  +  1) ; 
b2  =  valloc  (NULL,  p  +  D ; 
for  (1=0;  i<p+l;  i++) 
■{ 


(1  -  D); 


/*  clean  house  */ 

/*  Begin  Error  calculation  *//*R(l,:)  */ 
/*  R(p+1,:)  */ 


♦(bl  +  i)  =  R[0] [i];  /♦  R(l,:)  ♦/ 

b2[p  -  i]  =  R[p][i];  /♦  fliplr(  R(ptl,:)  )  ♦/ 


error  =  vdot  (bl,  ARparm,  p  +  1) 

vfree  (bl) ; 

vfree  (b2) ; 

inxfree  (X) ; 

inxfree  (R) ; 


vdot  (b2,  ARparm,  p  +  1); 


/♦  clean  house  ♦/ 


A.2.14  NLP.c 


#include  <stdlib.h> 

#include  "mathlib . h" 

#include  <math.h> 

#include  "slm_parameters_extern.h" 

Real_ Vector 

nip  (Real_Vector  r,  unsigned  int  r_length,  Real_Vector  mini,  Real_Vector  h, 
unsigned  int  h_length,  Real_Vector  bw,  Real_Vector  a, 
unsigned  int  a.length,  unsigned  int  g_length) 

■{ 


int  i,  j ,  k; 

Real.Natrlx  hargs, 
hlndx, 
tempmat ; 

Real_Vector  rho,  rho_l,  rho_2, 
hargCol , 
g> 

a.scaled; 

const  Real  inf  =  30000000, 

isinf  =  100000; 

Real  num; 


5  =  valloc  (NULL,  g_length) ;  /♦  allocate  vector  ♦/ 


/♦  arguments  for  the  "h"  function  ♦/ 

/♦  for  Indexes  of  h  function  ♦/ 

/♦  temp  matrix  ♦/ 

/♦  temp  usage  for  hargs  ♦/ 

/♦  hargs  column:  temporary  ♦/ 

/♦  vector  to  return  to  simulation. c  ♦/ 
/♦  Test  ♦/ 

Its  inf  in  matlab  ♦/ 

/♦  Min  value  for  "essentially  inf"  response  ♦/ 

/♦  temporary  storage  ♦/ 


rho  =  valloc  (NULL,  2  ♦  P  +  r_length) ; 
for  (i  =  0;  i  <  P;  1++) 


/♦  Concatenate  Initial  conditions  ft  rho  ♦/ 


rho[i]  =  0.0;  /♦  Initial  Values  of  zero  ♦/ 

rho[i  +  P  +  r_length]  =  inf;  /♦  Fixes  the  indexing  in  the  nonlinearity 

(forces  rhp(i>Nb)  =  Inf  (which  is  outside 
the  support  of  the  pdf  and  hence  hO  =  0! !  ♦/ 

} 

for  (i  =  P;  i  <  P  +  r_length;  i++) 
rho[i]  =  r[i  -  P] ; 

hargs  =  mxalloc  (NULL,  P  +  1,  Nb) ;  /♦  Initialize  the  args  for  "h"  function  ♦/ 

//tempmat  =  mxalloc(NULL,  P+1,  P+1);  //  allocated  within  Toeplitz.c 

rho_l  =  valloc  (NULL,  P  +  1) ; 

rho_2  =  valloc  (NULL,  P  +  1) ; 

a_scaled  =  valloc  (NULL,  a_length) ; 

a_scaled  =  vscale  (a,  a_length,  -1.0); 

for  (i  =  0;  i  <  Nb;  i++) 

■{ 


for  (j  =  0,  k  =  i  +  P;  j  <  P  +  1;  j++)  /♦  Allocate  rho(indx+P : indx+2P)  ♦/ 

■{  /♦  and  rho(indx+P: -1:  indx)  ♦/ 

rho_l[j]  =  rho[i  +  P  +  j]; 
rho_2[j]  =  rho[k  -  j] ; 

} 

tempmat  =  toeplitz  (rho_l,  rho_2,  P  +  1,  P  +  1);  /♦  Get  toeplitz  ♦/ 

hargCol  =  vmxmul  (tempmat,  a_scaled,  P  +  1,  P  +  1); 

for  (j  =  0;  j  <  P  +  1;  j++)  /♦  Copy  hargcol  into  appropriate  hargs  column  ♦/ 

■{ 


hargs[j][i]  =hargCol[j]; 

} 

mxf ree  (tempmat) ; 
vfree  (hargCol) ; 

/♦  Now  that  we  have  the  arguments  (a  matrix) ,  we  want  to  make  use  of  the 

histogram  estimate  of  the  pdf  (actually  the  derivative  of  the  natural  log  of 
the  pdf) .  So  we  take  the  "hargs"  matrix  and  change  it  to  INDEXES  of  the 
"h"  function,  utilizing  the  structure  of  the  pdf  (the  lowest  terms  are  on  the 
left  side  of  the  histogram). 

This  routine  uses  the  trick  of  assuming  the  minimum  values  in  the  histogram 
(mini)  is  located  in  the  first  bln,  and  then  we  use  the  combination  of  the 
breakpoint  information  and  the  floor  command  to  develop  an  Integer  that 
puts  the  value  of  "hargs"  in  the  correct  bin.  ♦/ 
tempmat  =  mxalloc  (NULL,  P  +  1,  Nb) ;  /♦  Allocate  matrix  of  ones  ♦/ 

mxinit  (tempmat,  P  +  1,  Nb,  1.0); 

hindx  =  mxsub  (hargs,  mxscale  (tempmat,  P  +  1,  Nb,  ♦mini),  P  +  1,  Nb) ; 
useinput_  =1;  /♦  overwrite  input  matrix  ♦/ 
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} 


hindx  =  mxscale  (hindx,  P  +  1,  Ifb,  (1.0  /  bw[0])); 
for  (1=0;  i  <  P  +  1;  i++)  /♦  floor(  hindx  )  */ 

for  (j  =  0;  j  <  Nh;  j++) 

hindxCi] [j]  =  floor  (hindxCi] [j] ) ; 
hindx  =  inxadd  (hindx,  tenpmat,  P  +  1,  Nh) ; 

useinput_  =0;  /♦do  not  overwrite  anymore  (DEFAULT) 

mxfree  (hargs);  //  printfC'l  "); 

mxfree  (tempmat);  //  printf ("2\n") ; 

/♦  Find  arguments  that  lie  outside  the  support  of  the  "h"  function  and  throw 
those  values  into  a  "junk"  bin  (assume  their  probability  Is  approx,  zero  ♦/ 
for  (1=0;  1  <  P  +  1;  i++)  /♦  Apply  "junk"  to  junkbin  ♦/ 

for  (j  =  0;  j  <  Nb;  j++) 

■{ 

if  ((hlndxCi]  [j]  <  1.0)  II  (hlndx[i][j]  >  (h_length  -  1))) 
hindxCi] [j]  =  (Real)  h_length; 

} 


/♦  ♦♦♦  Develop  non-linearity 

The  following  reshaping  makes  h(hlndx)  appear  In  the  correct  format 
We  want  P+1  rows  of  length  Nb,  where  the  arguments  for  "h"  In  each  row 
are  a  linear  combination  of  P+1  samples  of  the  rho  vector.  Each  row  Is  a 
shifted  version  to  accomplish  the  necessary  summations  for  g(rho)  ♦/ 
tempmat  =  mxalloc  (NULL,  P  +  1,  Nb) ;  /♦  h(hindx)  ♦/ 

for  (i  =  0;  i  <  P  +  1;  i++) 

■{ 

for  (j  =  0;  j  <  Nb;  j++) 

tempmat [1] [j]  =  h[(int)  (hindx[i][j]  -  1)]; 

} 


hargs  =  mxalloc  (NULL,  Nb,  P  +  1);  /♦  temp  matrix  for  reshaping  ♦/ 

hargs  =  mxtransp  (tempmat,  P  +  1,  Nb) ;  /♦  h(hindx’)  ♦/ 
mxfree  (hindx) ; 

hindx  =  mxalloc  (NULL,  P  +  1,  Nb) ; 
hindx  =  mxtransp  (hargs,  Nb,  P  +  1); 

/♦  Perform  a^reshape(h(hindx’) ,Nb,P+l) ’ ;  ♦/ 

for  (j  =  0;  j  <  Nb;  j++)  /♦  Multiply  over  each  column  ♦/ 

■{ 


} 


num  =0.0; 

for  (1=0;  i  <  P  +  1;  i++) 
num  =  num  +  (♦(a  +  i) 
*(g  +  j)  =  num; 


/♦  Initialize  num  to  zero  ♦/ 
/♦  Do  each  row  by  col  multiplication  ♦/ 
♦  hindxCi] [j] ) ; 


mxfree  (hindx) ; 
mxfree  (tempmat) ; 
mxfree  (hargs) ; 
vfree  (rho) ; 
vfree  (rho_l) ; 
vfree  (rho_2) ; 
return  g; 


♦/ 


A. 2. 15  Read  Farms. c 


/♦  This  function  reads  In  the  simulation  parameters  from  a  file 

♦  designated  on  the  command  line. 

♦ 

♦  The  file  should  have  one  parameter  per  line.  In  this  format; 

♦  <value>  <variable> 

♦  For  example,  512  G_b 

♦ 

♦  The  order  of  the  variables  should  be  as  follows; 

♦  ar_method(4) ,  iidtype(2),  Itype(2),  limit(4),  Bh,  Eb_max, 

♦  Eb_min,  I,  iidparms[2],  ISR,  it,  K,  L,  Nc,  P,  PG,  sigmat, 

♦  T,  z,  freq[2],  poly[2] 

*/ 

#include  <stdio.h> 

#include  <math.h> 

#include  <stdlib.h> 

#include  "mathlib . h" 

#include  "slm_parameters_extern.h" 

Intread.parms  (char  ♦filename) 

■{ 


int  i,  k, 

LI, 

numllne, 

space, 

char_variables  =  4, 
int_variables  =  16, 
numcommands  =  24, 
maxlength  =  20; 
char  commands [24] [20] ; 

float  j ; 

FILE  ♦parms; 

Real_Vector  poly_temp,  poly_temp2; 
Complex  initval  =  ■[ 

0,  0 

}; 


/♦  Counter  variables  ♦/ 

/♦  Lengths  for  recursive  poly[]  ♦/ 

/♦  Counter  for  current  line  being  processed  ♦/ 

/♦  Location  of  space  in  line  ♦/ 

/♦  number  of  variables  that  are  char  arrays  ♦/ 

/♦  number  of  variables  that  are  ints  ♦/ 

/♦  26  variable  assignments  ♦/ 

/♦  Max  length  of  a  line  ♦/ 

/♦  22  variable  assignments,  each 

♦  max  length  19  chars 

/♦  temp  storage  ♦/ 

/♦  For  opening  parameter  file  ♦/ 


/♦  For  vector  initialization  ♦/ 


if  ((parms  =  fopen  (filename,  "r"))  ==  NULL)  /♦  Open  parameter  file  ♦/ 
■[ 


1  for  NULL  ♦/ 
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printf  ("Caimot  open  parameter  file.Xn"); 
exit  (1) ; 

} 

for  (1=0;  1  <  numcommands;  i++)  /♦  Read  each  line  into  array  */ 

fgets  (commands [1] ,  (maxlength  -  1),  parms); 

/*  Input  character  arrays  first  */ 

for  (numline  =  0;  numline  <  char .variables ;  numline++) 

■{ 

space  =  0;  /*  Find  SPACE  character  */ 

for  (i  =  0;  i  <  maxlength;  i++) 

■{ 

if  (commands [numline] [1]  ==  ’  ’) 

■{ 

space  =  1; 


} 

} 

/*  Assign  string  to  defined  character  variables  */ 
switch  (numline) 

■{ 


} 


case 


case 


case 


case 


} 


0: 

for  (1=0;  i  <  space;  i++) 

ar .method [1]  =  commands [numline] [1] ; 

break; 

1: 

for  (1=0;  1  <  space;  i++) 

lldtype[l]  =  commands [numline] [1] ; 

break; 


for  (1=0;  i  <  space;  i++) 

Itype[l]  =  commands [numline] [1] ; 

break; 


for  (1=0;  i  <  space;  i++) 

limit [1]  =  commands [numline] [1] ; 

break; 


/*  Read  in  Integer  variables  */ 

for  (numline;  (numline  <  (char .variables  +  int.varlables)) ;  numline++) 

■[ 

i  =  atoi  (commands [numline]) ;  /*  Convert  char  number  to  int  */ 

/*  Assign  int  to  correct  variable  */ 
switch  (numline) 

■[ 

case  (4) ; 


Bh  =  1; 
break; 
case  (5) ; 

Eb.max  =  1; 
break; 
case  (6) ; 

Eb.min  =  1; 
break; 
case  (7) ; 

1  =  1; 
break; 
case  (8) ; 

iidparms[0]  =  1; 
break; 
case  (9) ; 

ildparms[l]  =  1; 
break; 
case  (lO) ; 

ISR  =  1; 
break; 
case  (ll) ; 

it  =  1; 
break; 
case  (12); 


K  =  1; 
break; 
case  (13) ; 

L  =  1; 
break; 
case  (14) ; 

Nc  =  1; 
break; 
case  (15) ; 

P  =  1; 
break; 
case  (16) ; 

PG  =  i; 
break; 
case  (17) ; 

slgmat  =  1; 
break; 
case  (18) ; 

T  =  1; 
break; 
case  (19) ; 
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z  =  i; 
break; 

} 

} 

/*  Read  in  float  variables  */ 

for  (nunline;  nunline  <  nuncomands;  nuiiillne++) 

■{ 

j  =  (float)  atof  (commands [numline] ) ;  /*  Convert  char  number  to  float  */ 

/*  Assign  float  to  correct  variable  */ 
switch  (numline) 

■{ 

case  (20) ; 

freq[o]  =  j;  /*  should  be  calculated,  not  input  */ 

freq[0]  =  z  *  Bh  *  Nc; 
break; 
case  (21) ; 

freq[l]  =  j;  /*  should  be  calculated,  not  input  */ 

freq[l]  =  z  *  Bh  *  Nc  *  0.3; 
break; 
case  (22) ; 

polynomial [o]  =  j; 

break;  /*  If  polynomial  length  >  2,  change  */ 

case  (23) :  /*  the  recursive  algorithm  below  */ 

polynomial [l]  =  j; 

break;  /*  as  needed.  */ 

} 

} 

/************************************** 

/*  Finish  Global  Assignments  */ 

/*  Spread  spectrum  parameters  */ 

Nb  =  PG  ♦  Nc; 

B  =  z  *  Bh; 

Ntot  =  B  *  Nb; 

/*  AR  Model  Parameters  */ 

Bar  =  Bh  /  64; 

Bg  =  Bh; 

/*  Signal  parameters  ft  Thermal  noise  */ 
slgmat2  =  pow  (sigmat,  2); 

Eb_sigma  =  valloc  (NULL,  (Eb_max  -  Eb_min  +1)); 
for  (i  =  0;  1  <  (Eb_max  -  Eb_mln  +  1);  i++) 

■{ 

j  =  (double)  (Eb_min  +  i)  /  10.0; 

Eb_slgma[i]  =  (Real)  pow  (10.0,  j); 

} 

/*  Nixed  jammer  */ 

ISRcw  =  ISRpb  =  ISR  /  2; 

/*  Create  poly  via  recursively  increasing  the  order  of  the  filter  to  2''L  */ 

Li  =  2;  /♦  Set  length  of  polynomial  vector  */ 

poly_temp2  =  valloc  (NULL,  (Li  +  Li  -  1)); 

poly_temp2[0]  =  polynomial [0] ; 

poly_temp2[l]  =  polynomial [1] ; 

for  (i  =  0;  1  <  L;  i++) 

■{ 

poly_temp  =  valloc  (NULL,  (Li  +  Li  -  D); 

conv  (poly_temp2,  Li,  poly_temp2,  Li,  poly_temp);  /*  poly_temp  =  conv(.)  ♦/ 

vfree  (poly_temp2) ; 

poly_temp2  =  valloc  (NULL,  Li  +  Li  -  1); 

vcopy  (poly_temp2,  poly_temp,  Li  +  Li  -  1);  /*  Copy  poly_temp  to  poly_temp2  ♦/ 

vfree  (poly_temp) ; 

Li  =  Li  +  Li  -  1;  /*  Set  to  new  length  */ 

} 

vfree  (poly_temp2) ; 

poly  =  valloc  (NULL,  Li);  /*  Create  global  vector  */ 

vcopy  (poly,  poly_temp2,  Li); 

a_t  =  valloc  (NULL,  Li) ; 

vcopy  (a_t,  poly,  Li); 

useinput_  =  1; 

vscale  (a_t,  Li,  -1.0); 

useinput_  =  0; 

if  (ltype[0]  ==  ’p’  ftft  Itype[l]  =  ’b’)  /*  Is  using  PB  interferer,  P  =  number  */ 

P  =  Li  -  1;  /*  of  coefficients  */ 

/*  End  Global  Assignments  */ 

/**************************************/ 
f close  (parms) ; 
return  1 ; 


A. 2. 16  Receiver. c 

#lnclude  <math.h> 

#include  <stdlib.h> 

#include  "mathlib . h" 

#include  "sim_parameters_extern.h" 
void 

receiver  (Real_Vector  signal,  Real_Vector  noise,  char  *limit, 

Real_Vector  rcvd_tr) 

■{ 

int  i;  /*  Counter  */ 
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for  (1=0;  1  <  Ntot;  i++) 

*(rcvd_tr  +  i)  =  *(slgnal  +  i)  +  *(noise  +  1); 
if  (limit [0]  ==  ’hO  /♦  If  limiter  is  a  HARD  LIMITER  ♦/ 

/*  aka  NOISE  BLANKER  */ 
for  (i  =  0;  1  <  Ntot;  i++) 

■{ 

if  ((♦(rcvd_tr  +  i)  >  T)  II  (♦(rcvd_tr  +  i)  <  -T)) 
♦(rcvd_tr  +  i)  =0; 

} 


/*  If  SOFT  LIMITER  */ 

for  (1=0;  1  <  Ntot;  i++) 

■{ 

if  ((♦(rcvd_tr  +1)  >  T)) 
♦(rcvd_tr  +  i)  =  T; 
else  if  (*(rcvd_tr  +  1)  <  -T) 
♦(rcvd_tr  +  i)  =  -T; 

} 


A. 2. 17  Reshape. c 

/♦♦♦This  Is  the  C  version  of  the  Matlab  ReshapeO  function^^^^ 

♦♦♦written  by  Tayo  Ihimoyan,  0l/2000^^^^^^/ 

#include  "cmt .h" 

#include  "mathlib . h" 

/♦ - ♦/ 

Real_Matrix 

reshape  (Real.Matrlx  M,  unsigned  ml,  unsigned  nl,  unsigned  m2,  unsigned  n2) 

■{ 

/♦ 


reshape ()  returns  the  m2-by-n2  matrix  M2  whose  elements 
are  taken  column-wise  from  M(ml-by-nl). 

On  exit,  the  output  matrix  M2[][]  is  the  value  of  this  routine. 

♦/ 

/♦ - 

unsigned  1,  j,  k,  1; 

Real_Matrix  M2; 
if  (M  ==  NULL) 

■{ 

matherr_  ("reshape",  E_NULLPTR) ; 
return  NULL; 

} 

/♦returns  an  error  message  If  m2^n2  is  not  equal  to  M  elements  ♦/ 
if  (ml  ♦  nl  !=  m2  ♦  n2) 

■{ 

matherr_  ("reshape",  E_NULLPTR) ; 
return  NULL; 

} 

if  ((M2  =  mxalloc  (NULL,  m2,  n2))  ==  NULL) 
return  NULL; 

k  =  0; 

1  =  0; 

for  (j  =  0;  j  <  nl;  j++) 

■{ 

for  (1=0;  1  <  ml;  i++) 

■{ 

if  (1  >=  m2) 

■{ 

1  =  0; 
k++; 

} 

M2[l][k]  =  M[i][j]; 

1++; 

} 

} 

return  M2; 

} 


A. 2. 18  Th  Inter. c 


/♦  Interferer  must  be  Initialized  to  zero  (or  contain  values  from  another 
Interference  type)  before  using  this  function 
♦/ 

#include  <stdlib.h> 

#include  <math.h> 

#include  <stdio.h> 

#include  <time.h> 

#include  "mathlib . h" 

#include  "sim_parameters_extern.h" 
voidth_inter  (Real_Vector  interferer) 

■{ 
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} 


int  k, 
nil, 

sig; 

Complex  pps, 
Iseq_conj ; 
double  whlte.nean, 
long  Itlne; 
Real_Vector  white, 


/♦  Counter  ♦/ 

/*  Mean  */ 

/*  Std  dev  ♦/ 

/*  Power  per  sample  In  the  sequence  */ 


white_varp,  white_vars,  c; 

/*  For  seeding  */ 

/*  White  noise  sequence  */ 
b;  /*  numerator  of  filter  */ 

Complex. Vector  Isequence; 

Real  normalize,  req;  /*  Required  power  level  */ 

if  (didtypeCo]  ==  ’gO  ftft  (iidtypeCl]  ==  ’s’))  /*  Gaussian  Case  */ 

■{ 


mu  =  lldparmsCo];  /*  Easier  reading  this  way  */ 

slg  =  lldparms[l]; 

white  =  valloc  (NULL,  Ntot) ; 

Itlme  =  time  (NULL) ; 

srand  ((unsigned)  Itlme  /  2); 

for  (k  =  0;  k  <  Ntot;  k++)  /*  Create  white  noise  sequence  */ 

♦(white  +  k)  =  slg  *  normal  (1,  0)  -  mu; 

whlte.mean  =  stats  (white,  Ntot,  ftwhite.varp,  ftwhite.vars) ;  /*  Get  actual  mean/var  */ 

for  (k  =  0;  k  <  Ntot;  k++)  /*  Make  white  EXACTLY  Gaussian  */ 

♦(white  +  k)  =  (♦(white  +  k)  -  white.mean)  /  sqrt  (white.varp) ; 

/♦  Build  the  Interference  signal  using  the  AR  model  as  a  filter  ♦/ 
b  =  valloc  (NULL,  l) ; 

c  =  192.6592;  /♦  Maximum  gain  in  poly  filter  (in  dB)  ♦/ 

b[o]  =  pow  (10. 0,  (-C  /  20.0));  /♦  Filter  with  no  zeros  ♦/ 
filter  (b,  1,  poly,  P  +  1,  white,  Ntot,  interferer) ; 


/♦  Find  the  total  power  per  sample  In  the  sequence,  and  normalize  the  Interferer  ♦/ 
Isequence  =  Cvalloc  (NULL,  Ntot) ; 

for  (k  =  0;  k  <  Ntot;  k++)  /♦  Copy  real  interferer  to  complex  Isequence  ♦/ 

■{ 


Isequence [k] .r  =  ♦(interferer  +  k); 

Isequence [k] . 1  =  0.0; 

} 

fft42  (Isequence,  Ntot,  -1);  /♦  Compute  fft  of  Isequence  ♦/ 

pps.i  =  0.0;  /♦  Initialize  power  per  sample  to  zero  ♦/ 

pps.r  =0.0; 

for  (k  =  0;  k  <  Ntot;  k++)  /♦  Compute  total  power  per  sample  ♦/ 

■{ 

Iseq.conj  =  Conjg  (♦(Isequence  +  k)); 

pps  =  Cadd  (pps,  Cmul  (Iseq.conj,  ♦(Isequence  +  k))); 

} 

pps.r  =  (Real)  pps.r  /  pow  (Ntot,  2); 

req  =  pow  (10,  (ISR  /  10));  /♦  Compute  required  power  level  ♦/ 

normalize  =  sqrt  ((Real)  req  /  pps.r);  /♦  Normalize  interferer  to  ISR  db  ♦/ 

useinput.  =1;  /♦  Perform  normalization  "in-place"  ♦/ 

interferer  =  vscale  (interferer,  Ntot,  normalize); 

useinput.  =0;  /♦  Reset  to  original  value  ♦/ 

vfree  (white) ; 
vfree  (b) ; 

Cvfree  (Isequence) ; 


A. 2. 19  Toeplitz.c 


/♦  Toeplitz 

♦  A  Toeplitz  matrix  Is  defined  by  one  column,  coll,  and  one  row,  rowl,  of  lengths  numrows 

♦  and  numcols,  respectively.  If  coll[0]  !=  rowl[0],  then  coll[0]  will  override. 

♦/ 

#include  "mathlib . h" 

#include  <stdlib.h> 

Real.Matrix 

toeplitz  (Real.Vector  coll,  Real.Vector  rowl,  unsigned  int  numrows, 
unsigned  Int  numcols) 


■{ 


register  int  column,  row;  /♦  counter  variables  ♦/ 

Real.Matrix  toepmat;  /♦  Output  matrix,  in  Toeplitz  form  ♦/ 

toepmat  =  mxalloc  (NULL,  numrows,  numcols);  /♦  Create  matrix  ♦/ 

for  (column  =  0;  column  <  numcols;  column++)  /♦  Fill  in  Toeplitz  form  into  matrix  ♦/ 
■{ 


if  (column  =  0) 
■{ 


/♦  Set  first  column  ♦/ 


} 


for  (row  =  0;  row  <  numrows;  row++) 

toepmat [row] [column]  =  ♦(coll  +  row);  /♦  Set  first  column  ♦/ 


} 


} 


} 


return  toepmat; 


toepmat [o] [column]  =  ♦(rowl  +  column);  /♦  Set  first  row  ♦/ 

for  (row  =  1;  row  <  numrows;  row++)  /♦  Use  shift  register  to  do  rest  of  rows  ♦/ 
toepmat [row] [column]  =  toepmat [row  -  1] [column  -  1]; 
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Appendix  B 


MATLAB®  Plotting  Script 


X  C  Slnulatlon  Plotting  Routine 
X 

X  Spetenber  16,  1999 
X  Michael  Banys 
Xc_strlng_creatlon; 

X  Load  strings  for  use  In  plots 

do.spectrum  =1;  X  Does  spectral  plot  If  do.spectrum  ==  1 
X  OPEN  FILES 

pbar  =  fopenC’pbar.txt’ , ’rO  5 
pblc  =  f openC’pblc .txt ’ , ’r’) ; 
pblo  =  fopenC’pblo.txt’ , ’rO ; 
pblw  =  fopenC’pblw.txt’ , ’r’) ; 
pbwc  =  f openC’pbwc  .txt ’ , ’rO  ; 
snr_vals  =  fopen(’snr_vals_.txt’ , ’r’) ; 

X  READ  DATA 

snr_llst  =  f scanf (snr_vals, ’Xf ’ , Inf ) ; 

Pb_ar_avg  =  fscanf (pbar, ’Xf’ , Inf ) ; 

Pb_lc_avg  =  fscanf (pblc, ’Xf’ , Inf ) ; 

Pb_lo_avg  =  fscanf (pblo, ’Xf’ , Inf ) ; 

Pb_lw_avg  =  fscanf (pblw, ’Xf’ , Inf ) ; 

Pb_wc_avg  =  fscanf (pbwc, ’Xf’ , Inf ) ; 

XXXXXXXXXXXXX  PBPLOT  xxxxxxxxxxxxxxxxx 

figure; 
hold  on; 

X  Nay  need  to  use  plot  &  loglO(.) 
title (strcat (plot_tltlel , INTERFERER,plot_tltle2) ) ; 
xlabeK’Slgnal-to-Themal  Noise  Ratio  (dB)0; 
ylabeK’Probablllty  of  Bit  Error  (10''x)O; 


plot(snr_llst,  loglO(Pb_ar_avg) ,  ’x-’) ; 
plot (snr_llst ,  loglO (Pb_wc_avg) ,  ’ : ’ ) ; 
plot(snr_llst,  loglO(Pb_lo_avg) ,  ’o’) j 
plot (snr_llst ,  loglO (Pb_lw_avg) ,  ’ — ’ ) ; 
plot (snr_llst ,  loglO (Pb_lc_avg) ,  ’ - ’ ) ; 


legend(gca,  ’ARLO’,  ’Correlator  on  whitened’,  ’LO  detector’,  ’LO  on  whitened’,  ’Linear  Correlator’); 
set(gca, ’YTlck’ , [-7  -6  -5  -4  -3  -2  -1  0  1]); 

message  =  sprlntf (’Make  note  of  all  the  Important  variables  (ISR,  Bh,  P,  etc.)’); 
dlsp (message) ; 
hold  off; 

X . 

If  do_spectrum  =  1 

XXXXXXXXXXXXXX  SPPLOT  xxxxxxxxxxxxxxxxxx 

lnterferer_out  =  f open(’ Interf erer_out .txt ’ , ’r’) ; 
nolse_out  =  f open(’nolse_out .txt ’ , ’r’) ; 
rcvd_tr_out  =  f open( ’rcvd_tr_out . txt’ , ’r ’) ; 
a_hat_out  =  f open(’a_hat_out .txt ’ , ’r’) ; 

Interferer  =  fscanf (lnterferer_out, ’Xf ’ , Inf ) ; 
noise  =  fscanf (nolse_out, ’Xf’ , Inf ) ; 
rcvd_tr  =  fscanf (rcvd_tr_ out , ’Xf’ , Inf ) ; 
a_hat  =f scanf (a_hat_out , ’Xf ’ , Inf ) ; 

Ptt  =  spectrum(nolse  -  Interferer) ;  X  Spectrum  of  thermal  noise 

Ptt  =  10*logl0(Ptt ( ; ,1)) ;  X  Remove  phase  Info  and  scale 

delta_f  =  l/(length(Ptt)  -  1);  X  For  appropriately  scaling  the  frequency  axis 

f  =  0: (length(Ptt)  -  1);  X  Obtain  frequency  components 

f  =  f  *  delta_f;  X  Scale  accordingly 

Pss  =  spectrum(nolse  -  rcvd_tr);  X  Spectrum  of  transmitted  signal 

Pss  =  10*logl0(Fss( ; ,1)) ;  X  Remove  phase  Info  and  scale 

Pll  =  spectrum(lnterferer) ;  X  Spectrum  of  Interference  signal 

Pll  =  10*logl0(Fll( ;  ,1)) ;  X  Remove  phase  Info  and  scale 

Prr  =  spectrum(rcvd_tr) ;  X  Spectrum  of  received  signal 

Prr  =  10*logl0(Frr( ; ,1)) ;  X  Remove  phase  Info  and  scale 

Pww  =  spectrum(fllter(-a_hat,  1,  rcvd_tr));  X  Spectrum  of  whitened  signal 
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Fww  =  10*logl0(Fww( ; ,1)) ;  X  Remove  phase  Info  and  scale 

figure; 

hold  on; 

title (strcat (plot_title3, INTERFERER,plot_title2) ) ; 
xlabeK’ Frequency  (with  respect  to  sampling  frequency)’); 
ylabeK’Fower  Magnitude  (dB)’); 
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